function dae = H2P_VPC(t,v,u,param,ctrl)
%input: t-time; v-states (diff(x)+alg(z)), u-inputs, param-constants
% algebraic equations are scalled by a factor of 1e-10
% antiwindup with backcalculation for econom, attemp and steam valve. 

nx = param(end);                               % number differential states
x  = v(1:nx);                                         % differential states
z  = v(nx+1:end);                                        % algebraic states
%% States
TE       = x(1);                              % [K], economizer temperature
TD       = x(2);                                     % [K] drum temperature
TSh      = x(3);                              % [K] superheater temperature
MD       = x(4);                                    % [kg] drum mass holdup
MSh      = x(5);         % [kg] superheatr holdup to calculate the pressure
w        = x(6);     % [Hz] frequency (grid). assume turbine speed is equal
TA       = x(7);                             % [K] attemperator temperature

% controller states. integrated error for CVs
interrMD  = x(8);                         % integrated error for drum level
interrTA  = x(9);       % integrated error for supeheated steam temperature
interrTG  = x(10);              % integrated error for flue gas temperature
interrPm  = x(11);                             % integrated error for power            
interryzV = x(12);        % integrated error for steam valve position (VPC)

%controller states. integrated error for MVs for antiwindup:backcalculation
interrmP  = x(13);                                  % inregrated error pump
interrmB  = x(14);                 % integrated error for economizer bypass
interrmA  = x(15);                      % integrated error for attemperator
interruzV = x(16);       % integrated error for steam valve form antiwindup

% algebraic states
TGE      = z(1);             % [K] flue gas at economizer outlet. algebraic
TGD      = z(2);                              % [K] flue gas at drum outlet
TGSh     = z(3);                                 % [K] flue gas at Shoutlet
pT       = z(4);                          % [bar] pressure at turbine inlet
Pm       = z(5);                            % [kW] power (negative of work) 

% fuel and valve inputs
mGdummy  = z(6);               % [kg/s] dummy state to output the fuel used
zVdummy  = z(7);             % [-] dummy state to output the valve position

%% Controller setpoint (new MVs), bias, gain, integral time. 
MDsp  = u(2);                                    % [kg] drum level setpoint
mP0   = ctrl(1);                                     %[kg/s] feedwater bias
KcD   = ctrl(2);                                        % proportional gain
tauID = ctrl(3);                                            % integral time

TAsp  = u(4);                 % [kg] superheated steam temperature setpoint
mA0   = ctrl(4);                                  %[kg/s] attemperator bias
KcA   = ctrl(5);                                        % proportional gain
tauIA = ctrl(6);                                            % integral time

TGEsp  = u(3);                         % [kg] flue gas temperature setpoint
mB0    = ctrl(7);                                %[kg/s] flue gas temp bias
KcTG   = ctrl(8);                                       % proportional gain
tauITG = ctrl(9);                                           % integral time

Pmsp     = u(1);                           % [kW] mechanical power setpoint
zV0      = ctrl(10);                          %[kg/s] flue gas - heat input
KczPm    = ctrl(11);                                    % proportional gain
tauIzPm  = ctrl(12);                                        % integral gain

zVsp   = u(5);                                    % Steam pressure setpoint
mG0    = ctrl(13);                                      % bias, steam valve
KcGz   = ctrl(14);                                      % proportional gain
tauIGz = ctrl(15);                                         % integral gain

tauTmB = tauITG;       % tracking time constant. eqaul to int time. for now
tauTmA = tauIA;        % tracking time constant. eqaul to int time. for now
tauTmP = tauID;        % tracking time constant. eqaul to int time. for now   
tauTzV = tauIzPm;      % tracking time constant. equal to int time. for now


%% Physical Constants 
Tref  = param(1);      % [K] reference temperature for enthalpy calculation
TrefB = param(2);              % [K] reference temperature for vaporization
dHvap = param(3);                 % [kJ/kg] vaporization heat at TD nominal
cpg   = param(4);                        % [kJ/kg/K] flue gas specific heat
cpW   = param(5);                           % [kJ/kg/K] water specific heat
cpSt  = param(6);                           % [kJ/kg/K] steam specific heat
R     = param(7);                       %[m3*bar/K/mol] ideal gas constant
Mwt   = param(8);                              % [kg/mol] water molar mass

%% Parameters. process variables kept constant. can be changed to states
T0    = param(9);             % [K] inlet temperature, equal to pump outlet
TG0   = param(10);                          % [K]flue gas inlet temperature
ME    = param(11);                                  % [kg] economizer holdup
pC    = param(12);                                %[bar] condenser pressure
%% Design parameters from solving steady-state model wrt nominal conditions
UAE   = param(13);                      % [kW/s] heat transf*area economizer
UAD   = param(14);                           % [kW/s] heat transf*area drum
UASh  = param(15);                    % [kW/s] heat transf*area superheater
CvD   = param(16);         % "valve coeff" for vapour mass flow out of drum
Kv    = param(17);                          % "valve coeff" for steam valve
VSh   = param(18);                                % [m3] superheater volume
phid  = param(19);             % mass flow coefficient at design conditions

%% Generator constants
Pm0   = param(20);                                         % [kW] base load
w0    = param(21);                                 % [Hz] nominal frequency
M0    = param(22);    %[pu MW/pu frequency/s] angular momentum-machine base
D0    = param(23);  %[pu MW/pu frequency/s] electrical damping-machine base
PL    = param(24);                                  % [kJ] load disturbance
%% Specific enthalpy
% water/steam
H0   = cpW*(T0-Tref);              % [kJ/kg], water inlet specific enthalpy
HE   = cpW*(TE-Tref);      % [kJ/kg], water specific enthalpy at economizer
HD   = cpW*(TrefB-Tref)+dHvap+cpSt*(TD-TrefB); % [kJ/kg] drum specific enth
HSh  = cpW*(TrefB-Tref)+dHvap+cpSt*(TSh-TrefB);  % [kJ/kg] Sh specific enth
HA   = cpW*(TrefB-Tref)+dHvap+cpSt*(TA-TrefB);  % [kJ/kg] att specific enth
% gas
HGD  = cpg*(TGD-Tref);%[kJ/kg], flue gas specific entlaphy at econom outlet
HGE  = cpg*(TGE-Tref);  %[kJ/kg], flue gas specific entlaphy at drum outlet
HGSh = cpg*(TGSh-Tref);   %[kJ/kg], flue gas specific entlaphy at SH outlet
HG0  = cpg*(TG0-Tref);    % [kJ/kg], flue gas specific enthalpy at SH inlet

%% 0. Gas side. MV for steam valve (VPC).
errP    = Pm-Pmsp;                          % deviation from power setpoint
zV      = zV0-KczPm*errP-KczPm/tauIzPm*interrPm-1/tauTzV*interruzV; 
zVp     = min(1,max(zV,0));                         % valve is opened [0,1]
erryzV  = zV-zVsp;                          % error from setpoint - for VPC 
erruzV  = zV-zVp;         % error from applied plant input - for antiwindup

mG      = mG0-KcGz*erryzV-KcGz/tauIGz*interryzV; %[kg/s] flue gass flowrate

%% 1. Pump  (steady-state)
% For now, initial conditions match pump outlet.
% we could introduce small holdups for the split and bypass
errM  = MD-MDsp;                           % error from drum level setpoint
errTA = TA-TAsp;                              % [K] error superheated steam
errTG = TGE-TGEsp;                               % [K] error flue gas steam

mP     = mP0-KcD*errM-KcD/tauID*interrMD-1/tauTmP*interrmP;     % pump flow
mPp    = max(0,mP);                                   % applied plant input
errmP  = mP-mPp;              % error of controller output from plant input

mA     = mA0-KcA*errTA-KcA/tauIA*interrTA-1/tauTmA*interrmA; %attemperator 
mAp    = max(0,mA);               % input applied to the plant. sat min = 0 
errmA  = mA-mAp;

mB     = mB0-KcTG*errTG-KcTG/tauITG*interrTG-1/tauTmB*interrmB; %bypass 
mBp    = max(0,mB);        % input applied to the plant. no upper sat limit
errmB  = mB-mBp;           % difference controller output and plant input

mE     = mPp-mAp-mBp; %[kg/s] mass flow rate at economizer inlet and outlet

%% 2. Economizer. Only energy balance and no mass balance(ct holdup)
DTE    = (TGD+TGE)/2-(TE+T0)/2;     % temperature difference:algebraic mean
dTEdt  = 1/ME*(mE*(T0-TE) + UAE*DTE/cpW);   % dynamic energy balance. steam
alghGE = (mG*(HGD-HGE))*1e-10 - (UAE*DTE)*1e-10;  % alg. energy balance gas

%% 3.Mixing after economizer
mMB    = mE+mBp;             % [kg/s] mass flow rate after economizer bypass
TMB    = (mE*TE+mBp*T0)/mMB;   % [K] temperature after mixing.  Not a state
HMB    = cpW*(TMB-Tref);                        % [kJ/kg] specific enthalpy
%aTMBdum = TMB*1e-10-TMBdum; % need scaling of TMB is an algebraic state
%% 4. Drum. Steam:Dynamic mass&energy balances. Gss: algebraic energy
A5=5.11564;   B5=1687.537;  C5=-42.98;           % Antoine parameters water
pD     = 10^(A5-B5/(TD+C5));                          % [bar] drum pressure 
pSh    = MSh/Mwt*R*TSh/VSh;                    % [bar] superheater pressure
dPDSh  = pD-pSh;                                % [bar] pressure difference
mD     = CvD*dPDSh;                               % [kg/s] flow out of drum
dMDdt  = mMB-mD;                               % mass balance. differential

DTD    = (TGSh+TGD)/2-(TMB+TD)/2;   % [K] temperature difference (alg mean)
dTDdt  = 1/(MD*cpW)*(-cpSt*TD*(mMB-mD)+mMB*HMB-mD*HD+UAD*DTD); % steam enrg
alghGD = (mG*(HGSh-HGD))*1e-10 - (UAD*DTD)*1e-10;      % gas energy balance

%% 5. Superheater.Steam:Dynamic mass&energy balances. Gss: algebraic energy

mT      = zVp*Kv*(pSh-pT);                   % [kg/s] steam valve mass flow
mSh     = mT-mAp;             % [kg/s] from attemperator static mass balance
dMShdt  = mD-mSh;            % mass balance superheater to get the pressure

DTSh    = (TG0+TGSh)/2-(TSh+TD)/2;      % [K] temperature difference (mean)
alghGSh = mG*(HG0-HGSh)*1e-10-UASh*DTSh*1e-10;    % energy balance gas side
dTShdt  = 1/(MSh*cpSt)*(-TSh*cpSt*(mD-mSh)+mD*HD-mSh*HSh+UASh*DTSh); %steam

%% 6. Attemperator. Static energy balance to calculate the temperature
%algTA   = (mSh*HSh+mA*H0)*1e-10-mT*HA*1e-10;
MA = 10; % [kg], attemperator holdup

dTAdt    = 1/(MA*cpSt)*(mSh*HSh+mAp*H0-mT*HA);

%% 7. Valve and turbine. Assume isentropic expansion
TT    = TA;                  % [K] isenthalpic valve. temp at turbine inlet
phi   = mT*sqrt(TT)/pT;             % mass flow coeff at current conditions
algT  = phi-phid;                  % turbine map. constant flow coefficient
TC    = TT*(pC/pT)^(R*1e5/(cpSt*Mwt*1e3)); %[K]condenser temperature if gas

W     = mT*cpSt*(TC-TT);                      % [kW] turbine work. negative
algPm = Pm+W;              % algebraic equation to get the mechanical power                                        % [kW] mechanical power

%% 8. Generator. Rotating masses inertia (M) and feedback from consumers(D)
% Assume generator frequency and machine frequency are equal. 
% D: Feedback from consumers. when the frequency drops, the consumers take 
% less power which decreases the load, so the frequency stabilizes at some
% new steady-state.
M     = M0*Pm0/w0;    % [kW/Hz]  angular momentum of the machine (absolute). For machine inertia
D     = D0*Pm0/w0;            % [kW/Hz] electrical damping coeff (absolute).
dwdt  = 1/M*(Pm-PL-D*(w-w0));            % frequency. differential equation

%% 9. Condenser&pump. to be added to close the cycle.

%% 10. dummy states for the controller outputs
mGalg = mG-mGdummy;
zValg = zVp-zVdummy;

%% System of differential and algebraic equations (DAE)
diff  = [dTEdt;dTDdt;dTShdt;dMDdt;dMShdt; dwdt;dTAdt];     % diff equations
alg   = [alghGE;alghGD;alghGSh; algT;algPm;mGalg;zValg];         % algebraic equations
erry  = [errM;errTA;errTG;errP;erryzV];%error from setpoint. make zero
erru  = [errmP;errmB; errmA; erruzV]; % error from applied input. make zero
dae   = [diff;erry;erru;alg];                             % function output

%  figure(1492)
%  plot(t,zV,'*')
%  hold on
 
%  figure(1504)
%  plot(t,mG,'*')
%  hold on

end


