AirTerminalDevice/model/DHU_3.py

import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt
from sklearn.metrics import r2_score,mean_absolute_error,mean_absolute_percentage_error
try:
    import plotnine as gg
except:
    pass

from .._model._base import BaseModel
from ..components.coil import CoolingCoil2,SteamCoilFs,SteamCoilFs2,SteamCoilFs3
from ..components.wheel import WheelS3
from ..components.mixed import Mixed

def record(name,var):
    pm.Deterministic(f'{name}_mu',var)

def observe(name,var,observed,sigma=1):
    mu    = pm.Deterministic(f'{name}_mu',var)
    sigma = pm.HalfNormal(f'{name}_sigma',sigma=sigma)
    pm.Normal(name,mu=mu,sigma=sigma,observed=observed)


class DHU_3(BaseModel):
    
    def __init__(self) -> None:
        super().__init__()
    
    def fit(
        self,
        input_data   : pd.DataFrame,
        observed_data: pd.DataFrame,
        rw_FA_val    : bool = False,
        plot_TVP     : bool = True
    ):
        with pm.Model() as self.MODEL_PYMC:
            param_prior = {
                'wheel_1'    : WheelS3.prior('wheel_1'),
                'wheel_2'    : WheelS3.prior('wheel_2'),
                'coil_2'     : CoolingCoil2.prior('coil_2'),
                'coil_3'     : CoolingCoil2.prior('coil_3'),
                'steamcoil_1': SteamCoilFs2.prior('steamcoil_1'),
                'steamcoil_2': SteamCoilFs.prior('steamcoil_2'),
                'F_air'      : {
                    'HzP_X' : pm.HalfNormal('F_air_HzP_X',sigma=1,initval=1),
                    'HzP_H' : pm.HalfNormal('F_air_HzP_H',sigma=1,initval=0.1),
                    'HzR_B' : pm.HalfNormal('F_air_HzR_B',sigma=1,initval=0.5),
                    'X_base': pm.TruncatedNormal('F_air_X_base',mu=0.5,sigma=0.2,lower=0,initval=0.5),
                    'H_base': pm.TruncatedNormal('F_air_H_base',mu=0.6,sigma=0.2,lower=0,upper=0.999,initval=0.6),
                    'B_base': pm.TruncatedNormal('F_air_B_base',mu=0.2,sigma=0.1,lower=0,initval=0.1),
                },
                'mixed_1':{},
                'mixed_2':{}
            }
            
            if rw_FA_val:
                N          = len(input_data)
                period     = 30
                n_segments = int(np.ceil(N/period))
                remainder  = N % period
                repeat     = [period] * (n_segments - 1) + ([remainder] if remainder != 0 else [])
                rw     = pm.GaussianRandomWalk(
                    'rw',sigma=0.1,init_dist=pm.Normal.dist(mu=0,sigma=0.3),shape=n_segments)
                param_prior['F_air']['val_rw']  = pm.Deterministic('F_air_val_rw',pt.repeat(rw,repeat))
                param_prior['F_air']['val_pct'] = pm.Beta('F_air_val_pct',alpha=8,beta=1,initval=0.9)
            
            res = DHU_3.model(
                Tin_F        = input_data.coil_1_ToutA.values,
                Hin_F        = input_data.coil_1_HoutA.values,
                HzP          = input_data.HzP.values,
                HzR          = input_data.HzR.values,
                coil_1_TinW  = input_data.coil_1_TinW.values,
                coil_2_TinW  = input_data.coil_2_TinW.values,
                coil_3_TinW  = input_data.coil_3_TinW.values,
                coil_1_Val   = input_data.coil_1_Val.values,
                coil_2_Val   = input_data.coil_2_Val.values,
                coil_3_Val   = input_data.coil_3_Val.values,
                wheel_1_TinR = input_data.wheel_1_TinR.values,
                wheel_2_TinR = input_data.wheel_2_TinR.values,
                mixed_1_TinM = input_data.mixed_1_TinM.values,
                mixed_2_TinM = input_data.mixed_2_TinM.values,
                mixed_1_HinM = input_data.mixed_1_HinM.values,
                mixed_2_HinM = input_data.mixed_2_HinM.values,
                engine       = 'pymc',
                param        = param_prior
            )
            observe('mixed_1_ToutA',res['mixed_1']['ToutA'],observed=observed_data.mixed_1_ToutA.values)
            observe('mixed_1_DoutA',res['mixed_1']['DoutA'],observed=observed_data.mixed_1_DoutA.values)
            observe('wheel_1_ToutC',res['wheel_1']['ToutC'],observed=observed_data.wheel_1_ToutC.values)
            observe('coil_2_ToutA',res['coil_2']['ToutA'],observed=observed_data.coil_2_ToutA.values)
            observe('coil_2_DoutA',res['coil_2']['DoutA'],observed=observed_data.coil_2_DoutA.values)
            observe('wheel_2_ToutP',res['wheel_2']['ToutP'],observed=observed_data.wheel_2_ToutP.values)
            observe('wheel_2_DoutP',res['wheel_2']['DoutP'],observed=observed_data.wheel_2_DoutP.values)
            observe('wheel_2_ToutR',res['wheel_2']['ToutR'],observed=observed_data.wheel_2_ToutR.values)
            observe('steamcoil_1_FP',res['steamcoil_1']['FP'],observed=observed_data.steamcoil_1_FP.values,sigma=1000)
            observe('steamcoil_1_Fs',res['steamcoil_1']['Fs'],observed=observed_data.steamcoil_1_Fs.values,sigma=20)
            observe('steamcoil_2_Fs',res['steamcoil_2']['Fs'],observed=observed_data.steamcoil_2_Fs.values,sigma=20)
            
            record('wheel_2_ToutC',res['wheel_2']['ToutC'])
            record('mixed_2_ToutA',res['mixed_2']['ToutA'])
            record('wheel_1_FaP',res['wheel_1']['FP'])
            record('wheel_1_FaR',res['wheel_1']['FR'])
            record('wheel_1_FaC',res['wheel_1']['FC'])
            record('mixed_1_FaA',res['mixed_1']['FA'])
            record('mixed_1_FaM',res['mixed_1']['FM'])
            record('F_air_S',res['Fa']['Fa_S'])
            record('F_air_H',res['Fa']['Fa_H'])
            record('F_air_X',res['Fa']['Fa_X'])
            
            self.param_posterior = pm.find_MAP(maxeval=50000,include_transformed=False)
            param_posterior_reorder = {'F_air':{}}
            for equp_name in param_prior.keys():
                param_posterior_reorder.setdefault(equp_name,{})
                for param_name,param_value in self.param_posterior.items():
                    if '__' in param_name:
                        continue
                    if param_name == 'F_air_val_rw':
                        param_value = np.median(param_value[-5:])
                    if param_name.startswith(equp_name):
                        param_name_adj = param_name.replace(f'{equp_name}_','')
                        param_posterior_reorder[equp_name][param_name_adj] = param_value
            self.record_model(
                model_name   = 'DHU',
                model        = param_posterior_reorder,
                train_data   = {'x':np.array([1])},
                train_metric = {'R2':1,'MAE':1,'MAPE':1}
            )
            
            # 样本内预测数据
            TVP_data = []
            for param_name in self.param_posterior.keys():
                if param_name.replace('_mu','') not in observed_data.columns:
                    continue
                TVP_data.append(
                    pd.DataFrame(
                        {
                            'param_name': param_name.replace('_mu',''),
                            'real'      : observed_data.loc[:,param_name.replace('_mu','')].values,
                            'pred'      : self.param_posterior[param_name]
                        }
                    )
                )
            self.TVP_data = pd.concat(TVP_data,axis=0)
            
            group_by_data = self.TVP_data.groupby(['param_name'])[['pred','real']]
            self.TVP_metric = (
                pd.concat(
                    [
                        group_by_data.apply(lambda dt:r2_score(dt.real,dt.pred)),
                        group_by_data.apply(lambda dt:mean_absolute_error(dt.real,dt.pred)),
                        group_by_data.apply(lambda dt:mean_absolute_percentage_error(dt.real,dt.pred)),
                    ],
                    axis=1
                )
                .set_axis(['R2','MAE','MAPE'],axis=1)
                .sort_values(by='R2',ascending=True)
            )
            if plot_TVP:
                gg.options.figure_size = (10,10)
                plot = (
                    self.TVP_data
                    .pipe(gg.ggplot)
                    + gg.aes(x='real',y='pred')
                    + gg.geom_point()
                    + gg.facet_wrap(facets='param_name',scales='free')
                    + gg.geom_abline(intercept=0,slope=1,color='red')
                )
                plot.show()
            
            return self
            
    def predict(self,input_data:pd.DataFrame) -> dict:
        param_posterior = self.model_info['model_DHU']
        res = DHU_3.model(
            Tin_F        = input_data.coil_1_ToutA.values,
            Hin_F        = input_data.coil_1_HoutA.values,
            HzP          = input_data.HzP.values,
            HzR          = input_data.HzR.values,
            coil_1_TinW  = input_data.coil_1_TinW.values,
            coil_2_TinW  = input_data.coil_2_TinW.values,
            coil_3_TinW  = input_data.coil_3_TinW.values,
            coil_1_Val   = input_data.coil_1_Val.values,
            coil_2_Val   = input_data.coil_2_Val.values,
            coil_3_Val   = input_data.coil_3_Val.values,
            wheel_1_TinR = input_data.wheel_1_TinR.values,
            wheel_2_TinR = input_data.wheel_2_TinR.values,
            mixed_1_TinM = input_data.mixed_1_TinM.values,
            mixed_2_TinM = input_data.mixed_2_TinM.values,
            mixed_1_HinM = input_data.mixed_1_HinM.values,
            mixed_2_HinM = input_data.mixed_2_HinM.values,
            engine       = 'numpy',
            param        = param_posterior
        )
        return res
    
    def predict_system(self,input_data:pd.DataFrame) -> pd.DataFrame:
        pred_res = self.predict(input_data)
        system_output = {}
        for equp_name,output_info in pred_res.items():
            for output_name,output_value in output_info.items():
                system_output[f'{equp_name}_{output_name}'] = output_value
        system_output       = pd.DataFrame(system_output)
        system_output['Fs'] = system_output.steamcoil_1_Fs + system_output.steamcoil_2_Fs
        return system_output
        
    
    @classmethod
    def model(
        cls,
        Tin_F,        # 前表冷后温度
        Hin_F,        # 前表冷后湿度
        HzP,          # 处理侧风机频率 
        HzR,          # 再生侧风机频率
        coil_1_TinW,  # 前表冷进水温度
        coil_2_TinW,  # 中表冷进水温度
        coil_3_TinW,  # 后表冷进水温度
        coil_1_Val,   # 前表冷阀门开度
        coil_2_Val,   # 中表冷阀门开度
        coil_3_Val,   # 后表冷阀门开度
        wheel_1_TinR, # 前转轮再生侧温度
        wheel_2_TinR, # 后转轮再生侧温度
        mixed_1_TinM, # 回风温度（处理侧）
        mixed_1_HinM, # 回风湿度（处理侧）
        mixed_2_TinM, # 补风温度（再生侧）
        mixed_2_HinM, # 补风湿度（再生侧）
        engine,
        param
    ) ->  dict:
        
        # 水的质量流量
        coil_2_FW = coil_2_Val / 100
        coil_3_FW = coil_3_Val / 100
        
        # 空气的质量流量
        F_air_HzP_H = param['F_air']['HzP_H']
        F_air_HzP_X = param['F_air']['HzP_X']
        F_air_HzP_S = F_air_HzP_H + F_air_HzP_X
        F_air_HzR_B = param['F_air']['HzR_B']
        
        F_air_S_base  = 1
        F_air_X_base  = param['F_air']['X_base']
        F_air_H_base  = param['F_air']['H_base']
        F_air_B_base  = param['F_air']['B_base']
        F_air_val_rw  = param['F_air'].get('val_rw',0)
        F_air_val_pct = param['F_air'].get('val_pct',0)
        
        F_air_X_base_adj = F_air_X_base + F_air_val_rw
        F_air_H_base_adj = F_air_H_base - F_air_val_rw * F_air_val_pct
        F_air_B_base_adj = F_air_B_base - F_air_val_rw * (1 - F_air_val_pct)
        
        Fa_S = F_air_S_base + F_air_HzP_S * (HzP / 50)            
        Fa_H = F_air_H_base_adj + F_air_HzP_H * (HzP / 50)            
        Fa_X = F_air_X_base_adj + F_air_HzP_X * (HzP / 50)
        Fa_B = F_air_B_base_adj + F_air_HzR_B * (HzR / 50)           
        Fa_P = Fa_B + Fa_X + Fa_H - Fa_S  
        
        wheel_1_FaP    = Fa_S - Fa_H
        wheel_1_FaC    = Fa_X - wheel_1_FaP
        wheel_1_FaR    = Fa_P
        wheel_2_FaP    = Fa_S
        wheel_2_FaC    = wheel_1_FaC
        wheel_2_FaR    = wheel_1_FaC
        mixed_1_FaM    = Fa_H
        mixed_1_FaA    = wheel_1_FaP
        mixed_2_FaM    = Fa_B
        mixed_2_FaA    = wheel_1_FaC
        coil_2_FaA     = Fa_S
        coil_3_FaA     = Fa_S
        steamcoil_1_Fa = Fa_P
        steamcoil_2_Fa = wheel_1_FaC
        
        # 前转轮
        wheel_1_res = WheelS3.model(
            TinP   = Tin_F,
            HinP   = Hin_F,
            FP     = wheel_1_FaP,
            TinR   = wheel_1_TinR,
            HinR   = 0,
            FR     = wheel_1_FaR,
            TinC   = Tin_F,
            HinC   = Hin_F,
            FC     = wheel_1_FaC,
            engine = engine,
            param  = param['wheel_1']
        )
        
        # 处理侧混风（回风）
        mixed_1_res = Mixed.model(
            TinA   = wheel_1_res['ToutP'],
            HinA   = wheel_1_res['HoutP'],
            FA     = mixed_1_FaA,
            TinM   = mixed_1_TinM,
            HinM   = mixed_1_HinM,
            FM     = mixed_1_FaM,
            engine = engine
        )
        
        # 中表冷
        coil_2_res = CoolingCoil2.model(
            TinA   = mixed_1_res['ToutA'],
            HinA   = mixed_1_res['HoutA'],
            FA     = coil_2_FaA,
            TinW   = coil_2_TinW,
            FW     = coil_2_FW,
            engine = engine,
            param  = param['coil_2']
        )
        
        # 后转轮
        wheel_2_res = WheelS3.model(
            TinP   = coil_2_res['ToutA'],
            HinP   = coil_2_res['HoutA'],
            FP     = wheel_2_FaP,
            TinC   = wheel_1_res['ToutC'],
            HinC   = wheel_1_res['HoutC'],
            FC     = wheel_2_FaC,
            TinR   = wheel_2_TinR,
            HinR   = 0,
            FR     = wheel_2_FaR,
            engine = engine,
            param  = param['wheel_2'],
        )
        
        # 后表冷
        coil_3_res = CoolingCoil2.model(
            TinA   = wheel_2_res['ToutP'],
            HinA   = wheel_2_res['HoutP'],
            FA     = coil_3_FaA,
            TinW   = coil_3_TinW,
            FW     = coil_3_FW,
            engine = engine,
            param  = param['coil_3']
        )
        
        # 后转轮湿度修正
        wheel_2_res_adj = WheelS3.model(
            TinP   = coil_2_res['ToutA'],
            HinP   = coil_2_res['HoutA'],
            FP     = wheel_2_FaP,
            TinC   = wheel_1_res['ToutC'],
            HinC   = wheel_1_res['HoutC'],
            FC     = wheel_2_FaC,
            TinR   = wheel_2_TinR,
            HinR   = wheel_2_res['HoutC'],
            FR     = wheel_2_FaR,
            engine = engine,
            param  = param['wheel_2'],
        )

        # 再生侧混风（排风）
        mixed_2_res = Mixed.model(
            TinA   = wheel_2_res_adj['ToutR'],
            HinA   = wheel_2_res_adj['HoutR'],
            FA     = mixed_2_FaA,
            TinM   = mixed_2_TinM,
            HinM   = mixed_2_HinM,
            FM     = mixed_2_FaM,
            engine = engine
        )
        
        # 前转轮湿度修正
        wheel_1_res_adj = WheelS3.model(
            TinP   = Tin_F,
            HinP   = Hin_F,
            FP     = wheel_1_FaP,
            TinR   = wheel_1_TinR,
            HinR   = mixed_2_res['HoutA'],
            FR     = wheel_1_FaR,
            TinC   = Tin_F,
            HinC   = Hin_F,
            FC     = wheel_1_FaC,
            engine = engine,
            param  = param['wheel_1']
        )
        
        # 前蒸气盘管
        steamcoil_1_res = SteamCoilFs2.model(
            TinA   = mixed_2_res['ToutA'],
            ToutA  = wheel_1_TinR,
            FA     = steamcoil_1_Fa,
            param  = param['steamcoil_1'],
            engine = engine
        )
        # steamcoil_1_res = SteamCoilFs3.model(
        #     TinA   = mixed_2_res['ToutA'],
        #     ToutA  = wheel_1_TinR,
        #     HinA   = mixed_2_res['DoutA'],
        #     HoutA  = mixed_2_res['HoutA'],
        #     FA     = steamcoil_1_Fa,
        #     param  = param['steamcoil_1'],
        #     engine = engine
        # )
        
        # 后蒸气盘管
        steamcoil_2_res = SteamCoilFs.model(
            TinA   = wheel_2_res_adj['ToutC'],
            ToutA  = wheel_2_TinR,
            FA     = steamcoil_2_Fa,
            param  = param['steamcoil_2'],
            engine = engine
        )
        
        return {
            'coil_2'     : coil_2_res,
            'coil_3'     : coil_3_res,
            'wheel_1'    : wheel_1_res_adj,
            'wheel_2'    : wheel_2_res_adj,
            'mixed_1'    : mixed_1_res,
            'mixed_2'    : mixed_2_res,
            'steamcoil_1': steamcoil_1_res,
            'steamcoil_2': steamcoil_2_res,
            'Fa':{
                'Fa_S': Fa_S,
                'Fa_H': Fa_H,
                'Fa_X': Fa_X,
                'Fa_P': Fa_P,
                'Fa_B': Fa_B,
            },
            'summary'    : {}
        }