function [dxdt]=xprime(t,x,flag,u,p) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [dxdt]=xprime_closed(t,x,flag,u,p)
% Reactor/separator process - closed loop simulations
% 
% NB:  Linearization
%
% Change: Controllers and setpoints  (Other controllers) 
%         Integral error (dIEdt)     
%         Stabilization control
%
% Written by: msg20000914 
%         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global Y tt YY type

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

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

ys=p(20:24); %Stabilization - setpoints

%States
Mri=x(1:2); Mri=Mri(:);
Mr=sum(Mri);
Mc1=x(3:2+NT);
Mc2=x(3+NT:2+2*NT);
Mc=Mc1+Mc2;
xcol=x(3:2*NT+2);
Mb=Mc(1); Md=Mc(NT);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Composition

%Recycle composition
xd1=Mc1(NT)/Mc(NT);
xd2=Mc2(NT)/Mc(NT);
xd=[xd1;xd2];

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

%Bottom product composition
xb1=Mc1(1)/Mc(1);
xb2=Mc2(1)/Mc(1);
xb=[xb1;xb2];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Control
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch type

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'stable'} %Conventional, self-optimizing, Luyben scheme

Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-30; tauI=1/6; Mrs=ys(3); IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr)+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'stable2'}  %Luyben Rule
   
Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
LT=LT+Kp*(Mds-Md); 

Kp=30; tauI=1/6; Mrs=ys(3); IE=x(NC+NC*NT+1);
R=R+Kp*(Mrs-Mr)+Kp/tauI*IE; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'stable3'}  %Luyben Scheme
   
Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
R=R+Kp*(Mds-Md); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'active'} 

Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-30; tauI=1/6; Mrs=ys(3); IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr)+Kp/tauI*IE;

Kp=12e4; tauI=2/3; xb2s=ys(4); IE=x(NC+NC*NT+2);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'L'} 

Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-10*10; tauI=1/6; Mrs=ys(3); %IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr);%+Kp/tauI*IE;

Kp=66.5845; tauI=0.1; xb2s=ys(4); IE=x(NC+NC*NT+1); 
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

LTs=ys(5);
LT=LTs;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'V'} 

Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-10*10; tauI=1/6; Mrs=ys(3); %IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr);%+Kp/tauI*IE;

Kp=-62.9624; tauI=0.1; xb2s=ys(4); IE=x(NC+NC*NT+1); 
LT=LT+Kp*(xb2s-xb(2))+Kp/tauI*IE;

VBs=ys(5);
VB=VBs;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'D'} 

Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
LT =LT+Kp*(Mds-Md);  

Kp=-30; tauI=1/6; Mrs=ys(3); IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr)+Kp/tauI*IE;

Kp=12e4; tauI=2/3; xb2s=ys(4); IE=x(NC+NC*NT+2);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Ds=ys(5);
D=Ds;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'F'}  %Luyben Rule

Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
LT=LT+Kp*(Mds-Md); 

Kp=30; tauI=1/6; Mrs=ys(3); IE=x(NC+NC*NT+1);  
R=R+Kp*(Mrs-Mr)+Kp/tauI*IE; 

Kp=12e4; tauI=2/3; xb2s=ys(4); IE=x(NC+NC*NT+2); 
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Fs=ys(5);
F=Fs;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'xd'} 

Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-30; tauI=1/6; Mrs=ys(3); IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr)+Kp/tauI*IE;

Kp=12e4; tauI=2/3; xb2s=ys(4); IE=x(NC+NC*NT+2);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Kp=7500; tauI=2/3; xd1s=ys(5); IE=x(NC+NC*NT+3);
LT=LT+Kp*(xd1s-xd(1))+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'xd2'} 

Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-10*10; tauI=1/6; Mrs=ys(3); %IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr);%+Kp/tauI*IE;

Kp=67.4862; tauI=0.1; xb2s=ys(4); IE=x(NC+NC*NT+1);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Kp=-51.4949; tauI=0.1; xd2s=ys(5); IE=x(NC+NC*NT+2);
LT=LT+Kp*(xd2s-xd(2))+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'L/F'} 

Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-30; tauI=1/6; Mrs=ys(3); IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr)+Kp/tauI*IE;

Kp=12e4; tauI=2/3; xb2s=ys(4); IE=x(NC+NC*NT+2);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

LTFs=ys(5);
LT=LTFs*F;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'V/F'} 

Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-10*10; tauI=1/6; Mrs=ys(3); %IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr);%+Kp/tauI*IE;

Kp=47.9983; tauI=0.2103; xb2s=ys(4); IE=x(NC+NC*NT+1);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Kp=67.2751; tauI=2.6046; VBFs=ys(5); IE=x(NC+NC*NT+2);
LT=LT+Kp*(VBFs-VB/F)+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'B/F'} 

Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
LT =LT+Kp*(Mds-Md);  

Kp=10; tauI=1/6; Mrs=ys(3); %IE=x(NC+NC*NT+1);
R=R+Kp*(Mrs-Mr);%+Kp/tauI*IE;

Kp=62.9625; tauI=0.1; xb2s=ys(4); IE=x(NC+NC*NT+1);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Kp=-1.0780; tauI=9.9922; BFs=ys(5); IE=x(NC+NC*NT+2);
F=F+Kp*(BFs-B/F)+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'D/F'} 

Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-10*10; tauI=1/6; Mrs=ys(3); %IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr);%+Kp/tauI*IE;

Kp=46.0137; tauI=4.0404; xb2s=ys(4); IE=x(NC+NC*NT+1);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Kp=-40.4118; tauI=5.0495; RFs=ys(5); IE=x(NC+NC*NT+2);
LT=LT+Kp*(RFs-R/F)+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'V/L'} 

Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-10*10; tauI=1/6; Mrs=ys(3); %IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr);%+Kp/tauI*IE;

Kp=83.1320; tauI=0.1; xb2s=ys(4); IE=x(NC+NC*NT+1);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Kp=-61.1160; tauI=0.1; VBLTs=ys(5); IE=x(NC+NC*NT+2);
LT=LT+Kp*(VBLTs-VB/LT)+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'B/L'} 

Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-10*10; tauI=1/6; Mrs=ys(3); %IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr);%+Kp/tauI*IE;

Kp=62.9626; tauI=0.1; xb2s=ys(4); IE=x(NC+NC*NT+1);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

BLTs=ys(5);
LT=B/BLTs;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'D/L'} 

Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-30; tauI=1/6; Mrs=ys(3); IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr)+Kp/tauI*IE;

Kp=12e4; tauI=2/3; xb2s=ys(4); IE=x(NC+NC*NT+2);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

RLTs=ys(5);
LT=R/RLTs;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'L/D'} 

Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-30; tauI=1/6; Mrs=ys(3); IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr)+Kp/tauI*IE;

Kp=12e4; tauI=2/3; xb2s=ys(4); IE=x(NC+NC*NT+2);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

LTRs=ys(5);
LT=LTRs*R;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'B/V'} 

Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-10*10; tauI=1/6; Mrs=ys(3); %IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr);%+Kp/tauI*IE;

Kp=0.4077; tauI=80.2227; xb2s=ys(4); IE=x(NC+NC*NT+1);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Kp=-0.1; tauI=61.6254; BVBs=ys(5); IE=x(NC+NC*NT+2);
LT=LT+Kp*(BVBs-B/VB)+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'D/V'} 

Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
LT=LT+Kp*(Mds-Md);  

Kp=-30; tauI=1/6; Mrs=ys(3); IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr)+Kp/tauI*IE;

Kp=12e4; tauI=2/3; xb2s=ys(4); IE=x(NC+NC*NT+2);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

RVBs=ys(5);
R=RVBs*VB;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'B/D'} 

Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=-10*10; tauI=1/6; Mrs=ys(3); %IE=x(NC+NC*NT+1);
F=F+Kp*(Mrs-Mr);%+Kp/tauI*IE;

Kp=78.8840; tauI=0.1000; xb2s=ys(4); IE=x(NC+NC*NT+1);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Kp=71.3969; tauI=0.1000; BRs=ys(5); IE=x(NC+NC*NT+2);
LT=LT+Kp*(BRs-B/R)+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'F/F0'} 

Kp=-10; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-10; Mds=ys(2);
LT =LT+Kp*(Mds-Md);  

Kp=10; tauI=1/6; Mrs=ys(3); %IE=x(NC+NC*NT+1);
R=R+Kp*(Mrs-Mr);%+Kp/tauI*IE;

Kp=62.9643; tauI=0.1; xb2s=ys(4); IE=x(NC+NC*NT+1);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Kp=1.3360; tauI=9.9660; FF0s=ys(5); IE=x(NC+NC*NT+2);
F=F+Kp*(FF0s-F/F0)+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'LS'}  %Luyben Scheme

Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Fs=ys(3);
F=Fs;

Kp=12e4; tauI=2/3; xb2s=ys(4); IE=x(NC+NC*NT+1);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Kp=7500; tauI=2/3; xd1s=ys(5); IE=x(NC+NC*NT+2);
LT=LT+Kp*(xd1s-xd(1))+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case{'BS'}  %Balanced Scheme

Kp=-30; Mbs=ys(1);
B=B+Kp*(Mbs-Mb);

Kp=-30; Mds=ys(2);
R =R+Kp*(Mds-Md);  

Kp=0.3e3; tauI=2/3; xrs=ys(3); IE=x(NC+NC*NT+1); %nominal 
F=F+Kp*(xrs-xr(1))+Kp/tauI*IE; 

Kp=12e4; tauI=2/3; xb2s=ys(4); IE=x(NC+NC*NT+2);
VB=VB+Kp*(xb2s-xb(2))+Kp/tauI*IE;

Kp=7500; tauI=2/3; xd1s=ys(5); IE=x(NC+NC*NT+3);
LT=LT+Kp*(xd1s-xd(1))+Kp/tauI*IE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
otherwise
   disp('No control'); 
end %switch 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%Column
D=R;ucol=[LT/60;VB/60;D/60;B/60;F/60;xr(1);1]; pcol=p(6:17);
xprime=colmmod(t,xcol,[],ucol,pcol); 

xprime=xprime*60;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Integral error
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch type
case {'unstable'}  
   dIEdt=[];                           %No integral error
case {'stable','stable2','stable3'}  
   dIEdt=[Mrs-Mr];     
case{'L','V','D','F','active'}
   dIEdt=[Mrs-Mr;xb2s-xb(2)];           %A,B,C,D
case{'xd'}
   dIEdt=[Mrs-Mr;xb2s-xb(2);xd1s-xd(1)];%E (conventional)
case{'xd2'}
   dIEdt=[xb2s-xb(2);xd2s-xd(2)];      %F
case{'L/F'}
   dIEdt=[Mrs-Mr;xb2s-xb(2)];          %G (self-optimizing)
case{'V/F'}
   dIEdt=[xb2s-xb(2);VBFs-VB/F];       %H
case{'B/F'}
   dIEdt=[xb2s-xb(2);BFs-B/F];         %I
case{'D/F'}
   dIEdt=[xb2s-xb(2);RFs-R/F];         %J
case{'V/L'}
   dIEdt=[xb2s-xb(2);VBLTs-VB/LT];     %K
case{'B/L'}
   dIEdt=[xb2s-xb(2);BLTs-B/LT];       %L
case{'D/L','L/D'}
   dIEdt=[Mrs-Mr;xb2s-xb(2)]; 
case{'B/V'}
    dIEdt=[Mrs-Mr;xb2s-xb(2)];        
case{'D/V'}
   dIEdt=[Mrs-Mr;xb2s-xb(2)];           %O
case{'B/D'}
   dIEdt=[xb2s-xb(2);BRs-B/R];         %P
case{'F/F0'}
   dIEdt=[xb2s-xb(2);FF0s-F/F0];       %Q
case{'LS'}
   dIEdt=[xb2s-xb(2);xd1s-xd(1)];      %R
case{'BS'}
   dIEdt=[xrs-xr(1);xb2s-xb(2);xd1s-xd(1)];      %xR
otherwise
   disp('Unknown alternative');
end %switch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%State variables
dxdt=[dMridt;xprime;dIEdt];
disp(sprintf('t=%0.4f',t));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Objective function
J=VB;
f=J;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Inequality constraints
g=zeros(6,1);
g(1)=Mr-2800;
g(2)=0.985-xb(2);
g(3)=-LT;
g(4)=-VB;
g(5)=-R;
g(6)=-B;

Ceq=[dxdt]; C=g;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Measurements
Y=[LT;VB;R;B;F;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
YY=[YY Y];
tt=[tt t];

