AirTerminalDevice/model/DHU/SDHU_AB.py

from typing import Union

import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt

from .._base._base_device import BaseDevice
from ...components import (
    coil_water,coil_steam,wheel2,wheel3,mixed
)
from ..utils.fit_utils import (
    observe,reorder_posterior
)
from ...tools.optimizer import optimizer
from ...tools.data_cleaner import DataCleaner



class SDHU_AB(BaseDevice):
    
    val_rw_adj_target = ('coil_2_ToutA','coil_2_DoutA')
    
    def __init__(
        self,
        DHU_type      = 'A',
        exist_Fa_H    = True,
        exist_Fa_H0   = False,
        wheel_1       = None,
        coolingcoil_2 = 'CoolingCoil2',
        heatingcoil_1 = 'SteamCoil',
        mixed_1       = 'Mixed',
        mixed_2       = 'Mixed',
        other_info    = None
    ) -> None:
        super().__init__()
        self.DHU_type = DHU_type.replace('SDHU_','')
        if self.DHU_type == 'A':
            wheel_1 = wheel_1 if wheel_1 is not None else 'WheelS3V3'
        elif self.DHU_type == 'B':
            wheel_1 = wheel_1 if wheel_1 is not None else 'WheelS2V2'
        else:
            raise Exception('SDHU_type must be A or B')
        
        self.components_str = {
            'wheel_1'      : wheel_1,
            'coil_2'       : coolingcoil_2,
            'heatingcoil_1': heatingcoil_1,
            'mixed_1'      : mixed_1,
            'mixed_2'      : mixed_2
        }
        self.exist_Fa_H = exist_Fa_H
        self.exist_Fa_H0 = exist_Fa_H0
        self.other_info = other_info if other_info is not None else {}
        self.record_load_info(
            components_str = self.components_str,
            DHU_type       = self.DHU_type,
            exist_Fa_H     = self.exist_Fa_H,
            exist_Fa_H0    = self.exist_Fa_H0,
            other_info     = self.other_info
        )
    
    @property
    def components(self):
        comp_map = {
            'WheelS2':wheel2,'WheelS3':wheel3,'CoolingCoil':coil_water,
            'SteamCoil':coil_steam,'Mixed':mixed
        }
        output ={}
        for comp_name,comp_model in self.components_str.items():
            if comp_model == 'SteamCoilVal':
                output[comp_name] = coil_steam.SteamCoilVal(
                    name     = comp_name,
                    Fs_rated = self.other_info[f'{comp_name}_Fs_rated']
                )
                continue
            for comp_map_k,comp_map_v in comp_map.items():
                if comp_model.startswith(comp_map_k):
                    output[comp_name] = getattr(comp_map_v,comp_model)(name = comp_name)
        return output
    
    @property
    def model_input_data_columns(self):
        columns = {
            'fan_1_Hz'    : 'fan_1_Hz',
            'fan_2_Hz'    : 'fan_2_Hz',
            'coil_2_TinW' : 'coil_2_TinW',
            'coil_2_Val'  : 'coil_2_Val',
            'wheel_1_TinR': 'wheel_1_TinR',
        }
        if self.exist_Fa_H:
            columns['mixed_1_TinM'] = 'mixed_1_TinM'
            columns['mixed_1_HinM'] = 'mixed_1_HinM'
        
        if self.exist_Fa_H0:
            columns['coil_1_ToutA'] = 'mixed_0_ToutA'
            columns['coil_1_HoutA'] = 'mixed_0_HoutA'
        else:
            columns['coil_1_ToutA'] = 'coil_1_ToutA'
            columns['coil_1_HoutA'] = 'coil_1_HoutA'
        
        if self.DHU_type == 'A':
            columns['Fa_TinA'] = 'Fa_TinA'
            columns['Fa_HinA'] = 'Fa_HinA'
        elif self.DHU_type == 'B':
            columns['mixed_2_TinM'] = 'mixed_2_TinM'
            columns['mixed_2_HinM'] = 'mixed_2_HinM'
            
        return columns
    
    @property
    def model_observe_data_columns(self):
        columns = {
            'mixed_1_ToutA': 'mixed_1_ToutA',
            'mixed_1_DoutA': 'mixed_1_DoutA',
            'coil_2_ToutA' : 'coil_2_ToutA',
            'coil_2_DoutA' : 'coil_2_DoutA',
        }
        if self.DHU_type == 'A':
            columns['wheel_1_ToutC'] = 'wheel_1_ToutC' # A类除湿机前转轮是三分转轮
        
        exclude_obs = self.other_info.get('exclude_obs',[])
        for col in exclude_obs:
            if col in columns:
                del columns[col]
        return columns
    
    def model(self,*args,**kwargs):
        if self.DHU_type == 'A':
            return model_A(*args,**kwargs)
        elif self.DHU_type == 'B':
            return model_B(*args,**kwargs)
        else:
            raise ValueError('DHU_type must be A or B')
    
    def fit(
        self,
        input_data   : pd.DataFrame,
        observed_data: pd.DataFrame,
        plot_TVP     : bool = True,
    ):
        if len(input_data) < 30:
            raise Exception('数据量过少')
        
        with pm.Model() as self.MODEL_PYMC:
            param_prior = {name:comp.prior() for name,comp in self.components.items()}
            if self.DHU_type == 'A':
                param_prior['F_air'] = AirFlow_SDHU_A.prior(exist_Fa_H = self.exist_Fa_H)
            elif self.DHU_type == 'B':
                param_prior['F_air'] = AirFlow_SDHU_B.prior(exist_Fa_H = self.exist_Fa_H)
            else:
                raise NotImplementedError
            
            res = self.model(
                **{k:input_data.loc[:,v].values for k,v in self.model_input_data_columns.items()},
                engine       = 'pymc',
                components   = self.components,
                param        = param_prior
            )
            for std_name,name in self.model_observe_data_columns.items():
                if name not in observed_data.columns:
                    raise Exception(f'Missing column: {name}')
                observed_data = observed_data.rename(columns={name:std_name})
                
                std_name_equp,std_name_point = std_name.rsplit('_',1)
                sigma = 1
                observe(
                    name     = std_name,
                    var      = res[std_name_equp][std_name_point],
                    observed = observed_data,
                    sigma    = sigma
                )
            
            self.param_posterior = pm.find_MAP(maxeval=50000,include_transformed=False)
        
        self.record_load_info(param_posterior = self.param_posterior)
        self.record_model(
            model_name   = 'ATD',
            model        = reorder_posterior(param_prior,self.param_posterior),
            train_data   = {
                'wheel_1_TinR': observed_data.loc[:,'wheel_1_TinR'].values,
                'fan_2_Hz'    : observed_data.loc[:,'fan_2_Hz'].values,
                'coil_2_DoutA': observed_data.loc[:,'coil_2_DoutA'].values,
            },
            train_metric = {'R2':1,'MAE':1,'MAPE':1}
        )
        self.TVP_data   = self.get_TVP(self.param_posterior,observed_data)
        self.TVP_metric = self.get_metric(self.TVP_data)
        if plot_TVP:
            self.plot_TVP(self.TVP_data).show()
        return self
    
    @property
    def F_air_val_rw(self):
        return None
    
    def set_F_air_val_rw(self,value:float):
        return self
    
    def clean_data(
        self,
        data         : pd.DataFrame,
        data_type    : list=['input','observed'],
        print_process: bool = True,
        fill_zero    : bool = False,
        save_log     : Union[str,None]  = None
    ) -> pd.DataFrame:
        data = data.replace(-9999,np.nan)
        clean_data = DataCleaner(data,print_process=print_process)
        
        filter_columns = []
        if 'input' in data_type:
            filter_columns += list(self.model_input_data_columns.values())
            clean_data = (
                clean_data
                .rm_rolling_fluct(window=60,fun='ptp',thre=0.1,include_cols=['State'])
                .rm_rule('State != 1')
                .rm_rule('fan_1_Hz < 10').rm_rule('fan_2_Hz < 10')
                .rm_outrange(method='raw',upper=140,lower=20,include_cols=['wheel_1_TinR'])
                .rm_outrange(method='raw',upper=15,lower=-60,include_cols=['coil_2_DoutA'])
            )
        if 'observed' in data_type:
            filter_columns += list(self.model_observe_data_columns.values())
        clean_data = clean_data.get_data(
            fill     = 0 if fill_zero else None,
            save_log = save_log
        )
        clean_data = clean_data.loc[:,filter_columns]
        return clean_data

    
    def optimize(
        self,
        cur_input_data: pd.DataFrame,
        wheel_1_TinR  : tuple = (70,120),
        fan_2_Hz      : tuple = (30,50),
        constrains    : list  = None,
        logging       : bool  = True,
        target        : str   = 'summary_Fs',
        target_min    : bool  = True
    ) -> list:
        constrains = [] if constrains is None else constrains
        cur_input_data = cur_input_data.iloc[[0],:]
        
        opt_var_boundary = {}
        if wheel_1_TinR is not None:
            opt_var_boundary['wheel_1_TinR'] = {'lb':min(wheel_1_TinR),'ub':max(wheel_1_TinR)}
        if fan_2_Hz is not None:
            opt_var_boundary['fan_2_Hz'] = {'lb':min(fan_2_Hz),'ub':max(fan_2_Hz)}
            
        opt_var_value = cur_input_data.loc[:,list(opt_var_boundary.keys())]
        oth_var_value = (
            cur_input_data
            .loc[:,list(set(self.model_input_data_columns.values()))]
            .drop(opt_var_value.columns,axis=1)
        )
        opt_res = optimizer(
            model            = self,
            opt_var_boundary = opt_var_boundary,
            opt_var_value    = opt_var_value,
            oth_var_value    = oth_var_value,
            target           = target,
            target_min       = target_min,
            constrains       = constrains,
            logging          = logging,
            other_kwargs     = {'NIND':2000,'MAXGEN':50}
        )
        return opt_res
    
    
    def plot_opt(
        self,
        cur_input_data: pd.DataFrame,
        target_min    : str   = 'summary_waste',
        coil_2_DoutA  : tuple = None
    ):
        if coil_2_DoutA is None:
            coil_2_DoutA = (
                self.model_info['model_train_info_ATD']['coil_2_DoutA_min'],
                self.model_info['model_train_info_ATD']['coil_2_DoutA_max']
            )
        data_input = (
            pd.MultiIndex.from_product(
                [
                    np.linspace(
                        self.model_info['model_train_info_ATD']['wheel_1_TinR_min']-5,
                        self.model_info['model_train_info_ATD']['wheel_1_TinR_max']+5,
                        1000
                    ),
                    np.linspace(
                        self.model_info['model_train_info_ATD']['fan_2_Hz_min']-5,
                        self.model_info['model_train_info_ATD']['fan_2_Hz_max']+5,
                        1000
                    ),
                ],
                names=['wheel_1_TinR','fan_2_Hz']
            )
            .to_frame(index=False)
        )
        for col in cur_input_data.columns:
            if col in data_input.columns:
                continue
            data_input[col] = cur_input_data.loc[:,col].iat[0]
        data_output = self.predict_system(data_input)
        data = (
            data_output
            .assign(
                wheel_1_TinR = data_input.loc[:,'wheel_1_TinR'],
                fan_2_Hz     = data_input.loc[:,'fan_2_Hz'],
            )
            .assign(coil_2_DoutA=lambda dt:dt.coil_2_DoutA.round(1))
            .loc[lambda dt:dt.coil_2_DoutA.between(*(min(coil_2_DoutA),max(coil_2_DoutA)))]
            .loc[lambda dt:dt.groupby('coil_2_DoutA')[target_min].idxmin()]
            .loc[lambda dt:dt.coil_2_DoutA.mod(0.5)==0]
        )
        import plotnine as gg
        plot = (
            data
            .pipe(gg.ggplot)
            + gg.aes(x='wheel_1_TinR',y='fan_2_Hz')
            + gg.geom_path(size=1)
            + gg.geom_point()
            + gg.geom_label(gg.aes(label='coil_2_DoutA'))
            + gg.geom_abline(slope=1,intercept=0,color='red',linetype='--')
        )
        return plot
    
    def plot_check(self,cur_input_data:pd.DataFrame) -> dict:
        pa1=self.curve(input_data=cur_input_data,x='wheel_1_TinR',y='wheel_1_DoutP')
        pa2=self.curve(input_data=cur_input_data,x='wheel_1_TinR',y='wheel_1_ToutP')
        pa3=self.curve(input_data=cur_input_data,x='wheel_1_TinR',y='wheel_1_EFF')
        pb1=self.curve(input_data=cur_input_data,x='fan_2_Hz',y='wheel_1_DoutP')
        pb2=self.curve(input_data=cur_input_data,x='fan_2_Hz',y='wheel_1_ToutP')
        pb3=self.curve(input_data=cur_input_data,x='fan_2_Hz',y='wheel_1_EFF')
        plot1 = (pa1|pa2|pa3)/(pb1|pb2|pb3)
        return {'plot1':plot1}
        
        
        
        
def model_A(
    Fa_TinA,        # 前表冷后温度
    Fa_HinA,        # 前表冷后湿度
    coil_1_ToutA,
    coil_1_HoutA,
    fan_1_Hz,     # 处理侧风机频率 
    fan_2_Hz,     # 再生侧风机频率
    coil_2_TinW,  # 中表冷进水温度
    coil_2_Val,   # 中表冷阀门开度
    wheel_1_TinR, # 前转轮再生侧温度
    engine    : str,
    components: dict,
    param     : dict,
    mixed_1_TinM = 0, # 回风温度（处理侧）
    mixed_1_HinM = 0, # 回风湿度（处理侧）
) ->  dict:
    
    # 水的质量流量
    coil_2_FW = coil_2_Val / 100
    # 空气的质量流量
    air_flow = AirFlow_SDHU_A.model(fan_1_Hz=fan_1_Hz,fan_2_Hz=fan_2_Hz,param=param)
    
    # 前转轮
    wheel_1_res = components['wheel_1'].model(
        TinP   = coil_1_ToutA,
        HinP   = coil_1_HoutA,
        FP     = air_flow['wheel_1_FaP'],
        TinR   = wheel_1_TinR,
        HinR   = 0,
        FR     = air_flow['wheel_1_FaR'],
        TinC   = Fa_TinA,
        HinC   = Fa_HinA,
        FC     = air_flow['wheel_1_FaC'],
        engine = engine,
        param  = param['wheel_1']
    )
    
    # 处理侧混风（回风）
    mixed_1_res = components['mixed_1'].model(
        TinA   = wheel_1_res['ToutP'],
        HinA   = wheel_1_res['HoutP'],
        FA     = air_flow['mixed_1_FaA'],
        TinM   = mixed_1_TinM,
        HinM   = mixed_1_HinM,
        FM     = air_flow['mixed_1_FaM'],
        engine = engine
    )
    
    # 中表冷
    coil_2_res = components['coil_2'].model(
        TinA   = mixed_1_res['ToutA'],
        HinA   = mixed_1_res['HoutA'],
        FA     = air_flow['coil_2_FaA'],
        TinW   = coil_2_TinW,
        FW     = coil_2_FW,
        engine = engine,
        param  = param['coil_2']
    )
    
    # 后转轮湿度修正
    wheel_1_res_adj = components['wheel_1'].model(
        TinP   = coil_1_ToutA,
        HinP   = coil_1_HoutA,
        FP     = air_flow['wheel_1_FaP'],
        TinR   = wheel_1_TinR,
        HinR   = wheel_1_res['HoutC'],
        FR     = air_flow['wheel_1_FaR'],
        TinC   = Fa_TinA,
        HinC   = Fa_HinA,
        FC     = air_flow['wheel_1_FaC'],
        engine = engine,
        param  = param['wheel_1']
    )
    
    # 前再生加热盘管
    heatingcoil_1_res = components['heatingcoil_1'].model(
        TinA   = wheel_1_res_adj['ToutC'],
        ToutA  = wheel_1_TinR,
        FA     = air_flow['heatingcoil_1_Fa'],
        param  = param['heatingcoil_1'],
        engine = engine
    )
    
    waste = cal_Q_waste(
        wheel_1_res       = wheel_1_res_adj,
        heatingcoil_1_res = heatingcoil_1_res,
        wheel_1_TinR      = wheel_1_TinR,
        fan_2_Hz          = fan_2_Hz
    )
    return {
        'coil_2'       : coil_2_res,
        'wheel_1'      : wheel_1_res_adj,
        'mixed_1'      : mixed_1_res,
        'heatingcoil_1': heatingcoil_1_res,
        'Fa'           : air_flow,
        'summary'      : {
            'Fs'  : heatingcoil_1_res['Fs'],
            **waste,
        }
    }
        
def model_B(
    mixed_2_TinM, # 再生侧入口温度
    mixed_2_HinM, # 再生侧入口湿度
    coil_1_ToutA,
    coil_1_HoutA,
    fan_1_Hz,     # 处理侧风机频率 
    fan_2_Hz,     # 再生侧风机频率
    coil_2_TinW,  # 中表冷进水温度
    coil_2_Val,   # 中表冷阀门开度
    wheel_1_TinR, # 前转轮再生侧温度
    engine    : str,
    components: dict,
    param     : dict,
    mixed_1_TinM = 0, # 回风温度（处理侧）
    mixed_1_HinM = 0, # 回风湿度（处理侧）
) ->  dict:
    # 水的质量流量
    coil_2_FW = coil_2_Val / 100
    # 空气的质量流量
    air_flow = AirFlow_SDHU_B.model(fan_1_Hz=fan_1_Hz,fan_2_Hz=fan_2_Hz,param=param)
    
    # 前转轮
    wheel_1_res = components['wheel_1'].model(
        TinP   = coil_1_ToutA,
        HinP   = coil_1_HoutA,
        FP     = air_flow['wheel_1_FaP'],
        TinR   = wheel_1_TinR,
        HinR   = mixed_2_HinM,
        FR     = air_flow['wheel_1_FaR'],
        engine = engine,
        param  = param['wheel_1']
    )
    # 处理侧混风（回风）
    mixed_1_res = components['mixed_1'].model(
        TinA   = wheel_1_res['ToutP'],
        HinA   = wheel_1_res['HoutP'],
        FA     = air_flow['mixed_1_FaA'],
        TinM   = mixed_1_TinM,
        HinM   = mixed_1_HinM,
        FM     = air_flow['mixed_1_FaM'],
        engine = engine
    )
    # 中表冷
    coil_2_res = components['coil_2'].model(
        TinA   = mixed_1_res['ToutA'],
        HinA   = mixed_1_res['HoutA'],
        FA     = air_flow['coil_2_FaA'],
        TinW   = coil_2_TinW,
        FW     = coil_2_FW,
        engine = engine,
        param  = param['coil_2']
    )
    # 前再生加热盘管
    heatingcoil_1_res = components['heatingcoil_1'].model(
        TinA   = mixed_2_TinM,
        ToutA  = wheel_1_TinR,
        FA     = air_flow['heatingcoil_1_Fa'],
        param  = param['heatingcoil_1'],
        engine = engine
    )
    waste = cal_Q_waste(
        wheel_1_res       = wheel_1_res,
        heatingcoil_1_res = heatingcoil_1_res,
        wheel_1_TinR      = wheel_1_TinR,
        fan_2_Hz          = fan_2_Hz
    )
    return {
        'coil_2'       : coil_2_res,
        'wheel_1'      : wheel_1_res,
        'mixed_1'      : mixed_1_res,
        'heatingcoil_1': heatingcoil_1_res,
        'Fa'           : air_flow,
        'summary'      : {
            'Fs'  : heatingcoil_1_res['Fs'],
            **waste,
        }
    }

class AirFlow_SDHU_A:
    
    @classmethod
    def model(cls,fan_1_Hz,fan_2_Hz,param):
        F_air_X2_base = 1     # 进入转轮处理侧的新风量
        F_air_H_base  = param['F_air'].get('H_base',0)
        F_air_P_base  = param['F_air']['P_base']
        
        Fa_H  = F_air_H_base + (fan_1_Hz/50) * param['F_air'].get('HzP_H',0)
        Fa_X2 = F_air_X2_base + (fan_1_Hz/50) * param['F_air']['HzP_X2']
        Fa_S  = Fa_H + Fa_X2
        Fa_P  = F_air_P_base + (fan_2_Hz/50) * param['F_air']['HzR_P']
        
        return {
            'wheel_1_FaP'     : Fa_X2,
            'wheel_1_FaC'     : Fa_P,
            'wheel_1_FaR'     : Fa_P,
            'coil_2_FaA'      : Fa_S,
            'mixed_1_FaM'     : Fa_H,
            'mixed_1_FaA'     : Fa_X2,
            'heatingcoil_1_Fa': Fa_P
        }
    
    @classmethod
    def prior(cls,exist_Fa_H):
        param = {}
        param['HzP_X2'] = pm.HalfNormal('F_air_HzP_X2',sigma=1,initval=0.5)
        param['HzR_P']  = pm.HalfNormal('F_air_HzR_P',sigma=1,initval=0.5)
        param['P_base'] = pm.TruncatedNormal('F_air_P_base',mu=1,sigma=0.2,lower=0,initval=1)
        if exist_Fa_H:
            param['H_base'] = pm.TruncatedNormal('F_air_H_base',mu=1,sigma=0.2,lower=0,initval=1)
            param['HzP_H']  = pm.HalfNormal('F_air_HzP_H',sigma=1,initval=0.5)
        return param
        
class AirFlow_SDHU_B:
    @classmethod
    def model(cls,fan_1_Hz,fan_2_Hz,param):
        F_air_X2_base = 1
        F_air_H_base  = param['F_air'].get('H_base',0)
        F_air_P_base  = param['F_air']['P_base']
        
        Fa_H  = F_air_H_base + (fan_1_Hz/50) * param['F_air'].get('HzP_H',0)
        Fa_X2 = F_air_X2_base + (fan_1_Hz/50) * param['F_air']['HzP_X2']
        Fa_S  = Fa_H + Fa_X2
        Fa_P  = F_air_P_base + (fan_2_Hz/50) * param['F_air']['HzR_P']
        return {
            'wheel_1_FaP'     : Fa_X2,
            'wheel_1_FaR'     : Fa_P,
            'coil_2_FaA'      : Fa_S,
            'mixed_1_FaM'     : Fa_H,
            'mixed_1_FaA'     : Fa_X2,
            'heatingcoil_1_Fa': Fa_P
        }
    @classmethod
    def prior(cls,exist_Fa_H):
        param = {}
        param['HzP_X2'] = pm.HalfNormal('F_air_HzP_X2',sigma=1,initval=0.5)
        param['HzR_P']  = pm.HalfNormal('F_air_HzR_P',sigma=1,initval=0.5)
        param['P_base'] = pm.TruncatedNormal('F_air_P_base',mu=1,sigma=0.2,lower=0,initval=1)
        if exist_Fa_H:
            param['H_base'] = pm.TruncatedNormal('F_air_H_base',mu=1,sigma=0.2,lower=0,initval=1)
            param['HzP_H']  = pm.HalfNormal('F_air_HzP_H',sigma=1,initval=0.5)
        return param
        

def cal_Q_waste(
    wheel_1_res,
    heatingcoil_1_res,
    wheel_1_TinR,
    fan_2_Hz
) -> dict:
    def waste_cond_func1(TinR):
        waste = 0.15 + 0.0001 * (TinR-70)**3
        return np.where(waste>0,waste,0)
    def waste_cond_func2(TinR):
        waste = 0.25 * (1 - np.exp(-0.04 * (TinR - 70)))
        return np.where(waste>0,waste,0)
    
    if isinstance(wheel_1_res['Qsen'],np.ndarray):
        WHERE = np.where
    else:
        WHERE = pm.math.switch
    
    heatingcoil_1_Q = heatingcoil_1_res['Q']
    heatingcoil_1_Q = WHERE(heatingcoil_1_Q>0,heatingcoil_1_Q,0)
    
    waste_Qsen1 = wheel_1_res['Qsen']
    waste_cond1 = heatingcoil_1_Q * waste_cond_func1(wheel_1_TinR)
    
    res = {
        'waste_Qsen1': waste_Qsen1,
        'waste_cond1': waste_cond1,
    }
    
    waste_out = heatingcoil_1_Q - wheel_1_res['Qsen'] - wheel_1_res['Qlat']
    waste_out = WHERE(waste_out>0,waste_out,0)
    waste = waste_Qsen1 + waste_cond1 + waste_out + fan_2_Hz/100
    res['waste_out']   = waste_out
    res['waste'] = waste
    
    return res