function [C, Ceq] = cola_lv_nonlcon(x,p,u0)

%******************************************************************
%                   GENERATING CONSTRAINTS

% Program for generating constraints for multi-effect distillation
% columns (consists of a HP and a LP column with integrated
% reboiler and condenser)

% Subroutine for the self optimisation program optmisation.m

% Calls model colamod1.m, the model for the HP column and colamod2.m
% the model for the LP column

%*******************************************************************
global P_HP P_LP T S

%x = S*x;         % Scaling back to "original" values

% Number of stages in columns:
NT = 16;                       % Stages in HP column
Nt = 26;                       % Stages in LP column
   
% Inputs, disturbances and states:
LT_HP = x(1);                  % Reflux in HP column
VB_HP = x(2);                  % Boilup in HP column
A     = x(3);                  % Area of reboiler/condenser (m^2)
X1 = x(4:3*NT+3);              % States HP column

LT_LP = x(3*NT+4);             % Reflux in LP column
VD_LP = x(3*NT+5);             % Condensed flow, LP column
X2 = x(3*NT+6:3*NT+5+3*Nt);    % States in LP column

zF_HP = p(1);                  % Light comp. fn. in feed HP col.
qF_HP = p(2);                  % Liq. fn. in feed to HP column
F_HP  = p(3);                  % Feed to HP column


% Max/min values for Inequality Constraints:
%--------------------------------------------

% HP Column
xD_HPmin = p(6);     % Min light component in HP column, p(6)
P_HPmax  = p(7);     % Max pressure HP column (Pa), p(7)
VB_HPmin = p(8);     % Min reboil flow (kmol/h), p(8) 
LT_HPmin = p(9);     % Min reflux flow (kmol/h), p(9)
VD_HPmin = p(10);    % Min. flow condensated (kmol/h) p(10)

% LP Column
xD_LPmin = p(11);    % Min light component in LP column, p(11)
P_LPmin  = p(12);    % Min pressure LP column (Pa), p(12)
VB_LPmin = p(13);    % Min reboil flow (kmol/h), p(13) 
LT_LPmin = p(14);    % Min reflux flow  (kmol/h), p(14)
VD_LPmin = p(15);    % Min. flow condensated (kmol/h) p(15)

% What about temp. constraints : Tc > Tr, Tsteam > Tr, Tc > Tcw

% P-controllers for Control of Reboiler and Condenser Hold-up:
%-------------------------------------------------------------
KcB_HP = -60; KcD_HP = -60;         % Controller gains (HP col)
KcB_LP = -60; KcD_LP = -60;         % Controller gains (HP col)

MDs_HP = 100; MBs_HP = 100;     % Nominal holdups (these are rather small)
MDs_LP = 100; MBs_LP = 100;     % Nominal holdups (these are rather small)
D0_HP = u0(1); B0_HP = u0(2);     % Nominal flows HP column (kmol/h)
D0_LP = u0(3); B0_LP = u0(4);     % Nominal flows LP column (kmol/h)

MB_HP = X1(NT+1);                    % Actual reboiler holdup (HP)
MD_HP = X1(2*NT);                    % Actual condenser holdup (HP)
D_HP = D0_HP + (MDs_HP-MD_HP)*KcD_HP;   % Distillate flow (HP)
B_HP = B0_HP + (MDs_HP-MB_HP)*KcD_HP;   % Bottoms flow (HP)

MB_LP = X2(Nt+1);                       % Actual reboiler holdup (LP)
MD_LP = X2(2*Nt);                       % Actual condenser holdup (LP)
D_LP = D0_LP + (MDs_LP-MD_LP)*KcD_LP;   % Distillate flow (LP)
B_LP = B0_LP + (MDs_LP-MB_LP)*KcD_LP;   % Bottoms flow (LP)

%---------------------------------------------------------------------

% Feed to LP column
zF_LP = X1(1);                 % Feed concn to LP col. (=xB_HP)
F_LP = B_HP;                   % Feedrate to LP column (kmol/h)
qF_LP =1;                      % Assume feed to LP is all liquid

% HP/LP Heat Exchanger Calculations:
%------------------------------------
Hvap_top_HP = 31710;            % From HYSYS at 5.2bar & 0.99 MeOH
Hvap_btm_LP = 40650;            % Assume pure water at 1 bar

Amax = 1000;                    % Max. area (m^2)
U=3000;                         % Overall h.t. coefficient (W/m^2 C)
DeltaT = X1(3*NT)-X2(2*Nt+1);   % (Tcond-Treb)
Q_HP = (U*A*DeltaT)*3600/1000;  % Duty in reboiler/condenser (kJ/h)

VD_HP = Q_HP/Hvap_top_HP;       % Condensed flow HP column (kmol/h)
VB_LP = Q_HP/Hvap_btm_LP;       % Boilup LP column (kmol/h)


%---------------------------------------------------------------

% Store all inputs and disturbances (goes into colamod.m)
U1(1)=LT_HP;U1(2)=VB_HP;U1(3)=D_HP;U1(4)=B_HP;U1(5)=F_HP;U1(6)=zF_HP;U1(7)=qF_HP;U1(8)=VD_HP;
U2(1)=LT_LP;U2(2)=VB_LP;U2(3)=D_LP;U2(4)=B_LP;U2(5)=F_LP;U2(6)=zF_LP;U2(7)=qF_LP;U2(8)=VD_LP;

% Call colamo1.m and colamod2.m
t = 0;
xprime1 = colamod1(t,X1,U1);
xprime2 = colamod2(t,X2,U2);

% Getting variables for contraints
% Pressure in columns
P_btm_HP = P_HP(1);      % Pressure in btm of HP column
P_top_HP = P_HP(NT);
P_btm_LP = P_LP(1);      % Pressure in btm of LP column
P_top_LP = P_LP(Nt);

% Fraction of MeOH in top of columns
xD_HP = X1(NT);     % Composition in top of HP column
xD_LP = X2(Nt);     % Composition in top of LP column
xB_HP = X1(1);      % Composition in btm of HP column
xB_LP = X2(1);      % Composition in btm of LP column

% Temperatures (K)
T_top_HP = X1(3*NT);        % Temperature in top of HP column
T_btm_LP = X2(2*Nt+1);      % Temperature in btm of LP column

T = [T_top_HP T_btm_LP];

% Generating Inequality (C<0) and Equality (Ceq) Constraints:
%------------------------------------------------------------
C_unscaled = [(P_top_HP-P_HPmax);(VB_HPmin-VB_HP); (VB_LPmin-VB_LP);... 
     (LT_HPmin-LT_HP); (LT_LPmin-LT_LP); (VD_HPmin-VD_HP);...
     (VD_LPmin-VD_LP)];    

% Scaling the inequality constraints
S = [ 1 1 1 1 1 1 1];
%S= [ 5e5 1000 1000 1000 1000 1000 1000];  %0.01 0.01 1e5 10
S1 = diag(S);

C = inv(S1)*C_unscaled;

% Scaling equality constraints
S2 = eye(NT*3+Nt*3+4)*100;
Ceq_unscaled = [xprime1; xprime2;(xD_HPmin-xD_HP); (xD_LPmin-xD_LP);(P_LPmin-P_top_LP);(A-Amax)];

Ceq = inv(S2)*Ceq_unscaled;
