function [C,Ceq,DC,DCeq]=nonlcon(ux,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [C,Ceq,DC,DCeq]=nonlcon3(ux,p)
% Reactor/separator process with recycle and inert
% Nominal optimization - constraints
%
% Change: Pure optimization, active constraint control, 
%         short-cut evaluation
%
% Written by :  msg20000914 
% Modified by:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global Y g %Ycol  V L x1 x2 x3 y1 y2 y3 g 

u=ux(1:6);
x=ux(7:size(ux,1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input
LT=u(1);
VB=u(2);
R =u(3);
B =u(4);
F =u(5);
F0=u(6);

%Parameter
VBmax=p(1);
x0=zeros(2,1); x0(1)=p(2); x0(2)=1-p(2);
k =p(4);
qv=p(5);
NC=p(6);
NT=p(7);
NF=p(8);
F0max=p(20);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%States

Mri=x(1:NC); Mri=Mri(:);
Mr=sum(Mri);
dMrdx=zeros(1,2*NT+8);
dMrdx(7:8)=1;

Mc1=x(3:NT+2);
Mc2=x(NT+3:2*NT+2);
Mc=Mc1+Mc2;
xcol=x(3:2*NT+2);

Mb=Mc(1);
dMbdx=zeros(1,2*NT+8); dMbdx([   9 NT+9])=1;
Md=Mc(NT);
dMddx=zeros(1,NT*NC+NC+6); dMddx([NT+8 2*NT+8])=1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Recycle composition
xd=[Mc1(NT)/Mc(NT);Mc2(NT)/Mc(NT)];
dxddx=zeros(2,2*NT+8); 
dxddx(1,NT+8)=1/Md;dxddx(2,2*NT+8)=1/Md;
dxddx(1,:)=dxddx(1,:)-Mc1(NT)/Md^2*dMddx;
dxddx(2,:)=dxddx(2,:)-Mc2(NT)/Md^2*dMddx;

%Product composition
xb=[Mc1(1)/Mc(1);Mc2(1)/Mc(1)];
dxbdx=zeros(2,2*NT+8); 
dxbdx(1,9)=1/Mb; dxbdx(2,NT+9)=1/Mb; 

dxbdx(1,:)=dxbdx(1,:)-Mc1(1)/Mb^2*dMbdx;
dxbdx(2,:)=dxbdx(2,:)-Mc2(1)/Mb^2*dMbdx;


%Reactor composition
xr(1:2)=Mri/Mr; xr=xr(:);

dxrdx=zeros(2,2*NT+8); dxrdx(:,7:8)=eye(2)/Mr; 
dxrdx(1,:)=dxrdx(1,:)-Mri(1)/Mr^2*dMrdx; 
dxrdx(2,:)=dxrdx(2,:)-Mri(2)/Mr^2*dMrdx;

%Reactor balances
v=[-1;1];
dMridt=F0*x0+R*xd-F*xr+Mr*xr(1)*k*v;


dfrdx=zeros(2,2*NT+8);
dfrdx(1,3)=xd(1); dfrdx(1,5)=-xr(1); dfrdx(1,6)=x0(1); dfrdx(1,7)=-k;
dfrdx(1,:)=dfrdx(1,:)+R*dxddx(1,:)-F*dxrdx(1,:);

dfrdx(2,3)=xd(2); dfrdx(2,5)=-xr(2); dfrdx(2,6)=x0(2); dfrdx(2,7)=k;
dfrdx(2,:)=dfrdx(2,:)+R*dxddx(2,:)-F*dxrdx(2,:);

dFdx=zeros(2*NT+2,2*NT+8);
dFdx(1:2,:)=dfrdx;

%Column
D=R;uxcol=[LT/60;VB/60;D/60;B/60;F/60;xr;1;xcol]; pcol=p(6:17);
[c,ceq,DC,DCeq]=nonlcon_colm(uxcol,pcol);xprime=ceq;

DCeq=DCeq'; 
DCeq(:,1:5)=DCeq(:,1:5)/60; %Time scale in column

dFdx(3:2*NT+2,[1:5])=DCeq(:,[1:5]);
%dFdx(3:2*NT+2,[9:2*NT+8])=DCeq(:,[10:2*NT+9]);
dFdx(3:2*NT+2,[9:2*NT+8])=DCeq(:,[9:2*NT+8]);
dFdx(3:2*NT+2,:)=dFdx(3:2*NT+2,:)+DCeq(:,6)*dxrdx(1,:);
dFdx(3:2*NT+2,:)=dFdx(3:2*NT+2,:)+DCeq(:,7)*dxrdx(2,:);

%df1=Mr-3000;      df1dx=dMrdx;
%df2=VB-1200;      df2dx=zeros(size(dMrdx)); df2dx(2)=1;
%df3=xb(2)-0.9895; df3dx=dxbdx(2,:);
df1=Mr-2800;      df1dx=dMrdx;
df2=VB-1400;      df2dx=zeros(size(dMrdx)); df2dx(2)=1;
df3=xb(2)-0.985;  df3dx=dxbdx(2,:);
 

%State variables
dxdt=[dMridt;ceq];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Objective function
J=-F0;%(0.1*B-0.05*F0-250*5/1e6*VB)*8760*50;
f=J;%/1e6;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Inequality constraints
g=zeros(10,1);
g(1)=Mr-2800;
g(2)=0.985-xb(2);
g(3)=-LT;
g(4)=-VB;
g(5)=-R;
g(6)=-B;
g(7)=-Mb;
g(8)=-Md;
g(9)= VB-VBmax;
g(10)=F0-F0max;

dg=zeros(10,2*NT+8);
dg(1,:)=dMrdx;
dg(2,:)=-dxbdx(2,:);
dg(3,1)=-1;
dg(4,2)=-1;
dg(5,3)=-1;
dg(6,4)=-1;
dg(7,[9 NT+9])     =-1;
dg(8,[NT+8 2*NT+8])=-1;
dg(9,2)=1;
dg(10,6)=1;


df1=Mr-2800;      df1dx=dMrdx;
df2=VB-1400;      df2dx=zeros(size(dMrdx)); df2dx(2)=1;
df3=xb(2)-0.985;  df3dx=dxbdx(2,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Pure optimization
Ceq=[dxdt]; C=[g];
DCeq=[dFdx'];  DC=[dg'];

%Include active constraints
%Ceq=[dxdt;df1;df3]; C=[g];
%DCeq=[dFdx' df1dx' df3dx'];  DC=[dg'];
%Ceq=[dxdt;df1;df2;df3]; C=[g];
%DCeq=[dFdx' df1dx' df2dx' df3dx'];  DC=[dg'];

%Short-cut evaluation of loss - few tight constraints
dysdx=zeros(size(df1dx));
%dys=[B-317.7517];  dysdx(4)=1;   %B
%dys=[VB-1.5747e+3];dysdx(2)=1; %VB
%dys=[LT-1.1278e+3];dysdx(1)=1; %LT
%dys=[R-414.6992];  dysdx(5)=1;   %R +
%dys=[P-32.2483]; dysdx(3)=1; %P 
%dys=[F-764.6992]; dysdx(6)=1; %F
%dys=[xr(1)-0.3091]; dysdx=dxrdx(1,:); %xr1 
%dys=[xr(2)-0.4264];dysdx=dxrdx(2,:); %xr2 
%dys=[xr(3)-0.2645]; dysdx=dxrdx(3,:); %xr3
%dys=[xb(1)-6.5148e-04]; dysdx=dxbdx(1,:); %xb1
%dys=[xb(3)-0.00984852342844]; dysdx=dxbdx(3,:); %xb3 
%dys=[xd(1)-0.5283];dysdx=dxddx(1,:); %xd1 
%dys=[xd(2)-0.0260]; dysdx=dxddx(2,:); %xd2
%dys=[xd(3)-0.4456]; dysdx=dxddx(3,:); %xd3
%dys=[R/F-414.6992/764.6992];dysdx(5)=1/F; dysdx(6)=-R/F^2; %R/F
%dys=[LT/F-1.1278e3/764.6992];dysdx(1)=1/F; dysdx(6)=-LT/F^2;%LT/F
%dys=[F/F0-764.6992/350];dysdx(6)=1/F0; %F/F0
%dys=[R/LT-414.6992/1.1278e3];dysdx(1)=-R/LT^2;dysdx(5)=1/LT; %R/LT
%Ceq=[dxdt;df1;df3;dys]; C=[g];
%DCeq=[dFdx' df1dx' df3dx' dysdx'];  DC=[dg'];

%C, max(C),pause
%Ceq, max(Ceq),pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Measurements
%Y=[LT;VB;R;B;F;Mr;Mb;Md;xr;xb;xd;J];%opt
Y=[LT;VB;R;B;F;F0;Mr;Mb;Md;xr;xb;xd;LT/F;VB/F;B/F;R/F;VB/LT;B/LT;R/LT;B/VB;R/VB;B/R;F/F0;J];%opt












