function sscon
%%% arrow.m needed
%%% http://www.mathworks.de/matlabcentral/fileexchange/278-arrow-m
% initial and final conditions :
x0 = [26; 16];        param.c10 = x0(1);  param.c20 = x0(2);   % [g/L]
xf = [130;6.8];       param.c1f = xf(1);  param.c2f = xf(2);   % [g/L]
param.V0 = 100;                                                      % [L]

% parameters : 
param.k       = 16.98;                                            % [m/h]                      
param.A       = 1;                                                % [m2]
param.clim    = 260;                                              % [g/L]
param.copt    = param.clim/exp(1);                                % [g/L]
param.K       = 1/(param.c10*param.V0);
param.P       = 0.95;
param.eps     = 1e-3;
param.N       = 100000;
param.wT      = 3.9;
param.wD      = 10;
param.mulcopt = param.clim*exp((param.wT - (param.wT^2 + 4*param.k*param.A*param.wT*param.wD)^(1/2))/(2*param.A*param.k*param.wD));

% control :
alpha0   = 0;
alpha1   = 1;
alphainf = 1000;

% Tolerances :
AbsTol = 1e-13;
RelTol = 1e-12;

% figures options : 
marker_size = 15;

% simulation step : 
step = 0.5;

%% minimum time : 

x01t      = [param.c10;param.c20;0];
options1t = odeset('Events',@event1t,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1t   = [0:step:param.N];
[time1t,state1t,te1t] = ode45(@model,tspan1t,x01t,options1t,alpha0,param);           

options2t = odeset('Events',@event2t,'AbsTol',AbsTol,'RelTol',RelTol);
x02t      = state1t(end,:);
tspan2t   = [time1t(end):step:param.N];
[time2t,state2t,te2t] = ode45(@model,tspan2t,x02t,options2t,alpha1,param);

options3t = odeset('Events',@event3t,'AbsTol',AbsTol,'RelTol',RelTol);
x03t      = state2t(end,:);
tspan3t   = [time2t(end):step:param.N];
[time3t,state3t,te3t] = ode45(@model,tspan3t,x03t,options3t,alpha0,param);

%% minimum diluant :

x01d      = [param.c10;param.c20;0];
options1d = odeset('Events',@event1d,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1d   = [0:step:param.N];
[time1d,state1d,te1d] = ode45(@model,tspan1d,x01d,options1d,alpha0,param);           

options2d = odeset('Events',@event2d,'AbsTol',AbsTol,'RelTol',RelTol);
x02d      = state1d(end,:);
tspan2d   = [time1d(end):step:param.N];
[time2d,state2d,te2d] = ode45(@model,tspan2d,x02d,options2d,alpha1,param);

options3d = odeset('Events',@event3d,'AbsTol',AbsTol,'RelTol',RelTol);
x03d      = state2d(end,:);
tspan3d   = [time2d(end):step:param.N];
[time3d,state3d,te3d] = ode45(@model,tspan3d,x03d,options3d,alphainf,param);

%% real minimum diluant :
wTr = 0.01;      % euro/h 
wDr = 1;         % euro/m3 
param.coptrmd = param.clim*exp((wTr - (wTr^2 + 4*param.k*param.A*wTr*wDr)^(1/2))/(2*param.A*param.k*wDr));

x01r      = [param.c10;param.c20;0];
options1r = odeset('Events',@event1r,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1r   = [0:step:param.N];
[time1r,state1r,te1r] = ode45(@model,tspan1r,x01r,options1r,alpha0,param);           

options2r = odeset('Events',@event2r,'AbsTol',AbsTol,'RelTol',RelTol);
x02r      = state1r(end,:);
tspan2r   = [time1r(end):step:param.N];
[time2r,state2r,te2r] = ode45(@model,tspan2r,x02r,options2r,alpha1,param);

options3r = odeset('Events',@event3r,'AbsTol',AbsTol,'RelTol',RelTol);
x03r      = state2r(end,:);
tspan3r   = [time2r(end):step:param.N];
[time3r,state3r,te3r] = ode45(@model,tspan3r,x03r,options3r,alphainf,param);

%% multi-objective : 

x01m      = [param.c10;param.c20;0];
options1m = odeset('Events',@event1m,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1m   = [0:step:param.N];
[time1m,state1m,te1m] = ode45(@model,tspan1m,x01m,options1m,alpha0,param);           

options2m = odeset('Events',@event2m,'AbsTol',AbsTol,'RelTol',RelTol);
x02m      = state1m(end,:);
tspan2m   = [time1m(end):step:param.N];
[time2m,state2m,te2m] = ode45(@model,tspan2m,x02m,options2m,alpha1,param);

options3m = odeset('Events',@event3m,'AbsTol',AbsTol,'RelTol',RelTol);
x03m      = state2m(end,:);
tspan3m   = [time2m(end):step:param.N];
[time3m,state3m,te3m] = ode45(@model,tspan3m,x03m,options3m,alphainf,param);

%% almost minimum time :

wTamt = 1;      % euro/h 
wDamt = 0.1;    % euro/m3 
param.coptamt = param.clim*exp((wTamt - (wTamt^2 + 4*param.k*param.A*wTamt*wDamt)^(1/2))/(2*param.A*param.k*wDamt));

x01amt      = [param.c10;param.c20;0];
options1amt = odeset('Events',@event1amt,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1amt   = [0:step:param.N];
[time1amt,state1amt,te1amt] = ode45(@model,tspan1amt,x01amt,options1amt,alpha0,param);           

options2amt = odeset('Events',@event2amt,'AbsTol',AbsTol,'RelTol',RelTol);
x02amt      = state1amt(end,:);
tspan2amt   = [time1amt(end):step:param.N];
[time2amt,state2amt,te2amt] = ode45(@model,tspan2amt,x02amt,options2amt,alpha1,param);

options3amt = odeset('Events',@event3amt,'AbsTol',AbsTol,'RelTol',RelTol);
x03amt      = state2amt(end,:);
tspan3amt   = [time2amt(end):step:param.N];
[time3amt,state3amt,te3amt] = ode45(@model,tspan3amt,x03amt,options3amt,alphainf,param);

%% figures :
marker_sizes = 4;
figure
plot([state1d(:,1);state2d(:,1);state3d(:,1)],[state1d(:,2);state2d(:,2);state3d(:,2)],'r-','MarkerEdgeColor','r','MarkerFaceColor','r','markersize',marker_sizes)               % minimum diluant
hold on
plot([state1r(:,1);state2r(:,1);state3r(:,1)],[state1r(:,2);state2r(:,2);state3r(:,2)],'co-','MarkerEdgeColor','c','MarkerFaceColor','c','markersize',marker_sizes)              % realistic minimum diluant
plot([state1m(:,1);state2m(:,1);state3m(:,1)],[state1m(:,2);state2m(:,2);state3m(:,2)],'kd-','MarkerEdgeColor','k','MarkerFaceColor','k','markersize',marker_sizes)              % economically optimal
plot([state1amt(:,1);state2amt(:,1);state3amt(:,1)],[state1amt(:,2);state2amt(:,2);state3amt(:,2)],'gx-','MarkerEdgeColor','g','MarkerFaceColor','g','markersize',marker_sizes)  % almost minimum time
plot([state1t(:,1);state2t(:,1);state3t(:,1)],[state1t(:,2);state2t(:,2);state3t(:,2)],'bs-','MarkerEdgeColor','b','MarkerFaceColor','b','markersize',marker_sizes)              % minimum time
plot(param.c10,param.c20,'go','MarkerEdgeColor','g','markersize',marker_size)
plot(param.c1f,param.c2f,'rx','markersize',marker_size,'MarkerEdgeColor','r','MarkerFaceColor','r')

vs = 18;
arrow([state2d(end,1);state2d(end,2)],[state3d(end,1);state3d(end,2)],vs,'Ends','stop','FaceColor','r','EdgeColor','r','LineStyle','--')
arrow([state2r(end,1);state2r(end,2)],[state3r(end,1);state3r(end,2)],vs,'Ends','stop','FaceColor','c','EdgeColor','c','LineStyle','--')
arrow([state2m(end,1);state2m(end,2)],[state3m(end,1);state3m(end,2)],vs,'Ends','stop','FaceColor','k','EdgeColor','k','LineStyle','--')
arrow([state2amt(end,1);state2amt(end,2)],[state3amt(end,1);state3amt(end,2)],vs,'Ends','stop','FaceColor','g','EdgeColor','g','LineStyle','--')
arrow([state2t(end,1);state2t(end,2)],[state3t(end,1);state3t(end,2)],vs,'Ends','stop','FaceColor','b','EdgeColor','b')
xlabel('c1'); ylabel('c2')
legend('minimum diluant','realistic minimum diluant','economically optimal','almost minimum time','minimum time','Location','SouthEast')
axis([0 300 4 17])
set(gcf, 'PaperUnits', 'centimeters')
set(gcf, 'PaperPosition', [3 10 16 10])
print -depsc 'soybean_state_fixed.eps'

figure
boundup = 1.1;
% minimum diluant : 
a1d = alpha0*ones(length(time1d),1);
a2d = alpha1*ones(length(time2d(1:3)),1);
a2d2= alpha1*ones(length(time2d(5:7)),1);
a2d3= alpha1*ones(length(time2d(3:5)),1);
plot([time1d;time2d(1:3)],[a1d;a2d],'r','MarkerEdgeColor','r','MarkerFaceColor','r','markersize',marker_sizes)
hold on   
% realistic minimum diluant
a1r = alpha0*ones(length(time1r),1);
a2r = alpha1*ones(length(time2r),1);
a3r = alphainf*ones(length(time3r),1);
plot([time1r;time2r;time3r],[a1r;a2r;a3r],'co-','MarkerEdgeColor','c','MarkerFaceColor','c','markersize',marker_sizes)
% economically optimal : 
a1m = alpha0*ones(length(time1m),1);
a2m = alpha1*ones(length(time2m),1);
a3m = alphainf*ones(length(time3m),1);
plot([time1m;time2m;time3m],[a1m;a2m;a3m],'kd-','MarkerEdgeColor','k','MarkerFaceColor','k','markersize',marker_sizes)
%almost minimum time :
a1amt = alpha0*ones(length(time1amt),1);
a2amt = alpha1*ones(length(time2amt),1);
a3amt = alphainf*ones(length(time3amt),1);
plot([time1amt;time2amt;time3amt],[a1amt;a2amt;a3amt],'gx-','MarkerEdgeColor','g','MarkerFaceColor','g','markersize',marker_sizes)
% minimum time : 
a1t = alpha0*ones(length(time1t),1);
a2t = alpha1*ones(length(time2t),1);
a3t = alpha0*ones(length(time3t),1);
plot([time1t;time2t;time3t],[a1t; a2t; a3t],'bs-','MarkerEdgeColor','b','MarkerFaceColor','b','markersize',marker_sizes)

plot([13.63;13.63],[1;1.1],'r',13.63,1.1,'^','MarkerEdgeColor','r','MarkerFaceColor','r')
plot(time2d(5:7),a2d2,'r','MarkerEdgeColor','r','MarkerFaceColor','r','markersize',marker_sizes)
plot(time2d(3:5),a2d3,'k.')
plot(time2m(end),boundup,'^','MarkerEdgeColor','k','MarkerFaceColor','k')
plot(time2amt(end),boundup,'^','MarkerEdgeColor','g','MarkerFaceColor','g')
plot(time2r(end),boundup,'^','MarkerEdgeColor','c','MarkerFaceColor','c')
arrow([7;0.2],[4.3854;0],1,'Ends','stop','FaceColor','k','EdgeColor','k','LineStyle','-','TipAngle',0.003,'Length',12)
text(7.2, 0.23  , 'mintf', 'Color', 'k');
xlabel('time'); ylabel('alpha')
legend('minimum diluant','realistic minimum diluant','economically optimal','almost minimum time','minimum time','Location','East')
axis([0 14 -0.1 1.12])
set(gcf, 'PaperUnits', 'centimeters')
set(gcf, 'PaperPosition', [3 10 16 10])
print -depsc 'soybean_control_fixed.eps'

function fun = model(t,x,u,param)
  c1 = x(1); c2 = x(2); alpha = u;
  q = param.k*param.A*log(param.clim/c1);
  
  dc1dt = param.K*c1^2*q*(1 - alpha);
  dc2dt = -c1*c2*param.K*q*alpha;
  daqdt = alpha*q;
  fun = [dc1dt;dc2dt;daqdt];

function [value,isterminal,direction] = event1amt(t,y,u,param)
  value      = y(1) - param.coptamt;
  isterminal = 1; direction  = 1;
function [value,isterminal,direction] = event1d(t,y,u,param)
  value      = y(1) - (param.clim - param.eps);
  isterminal = 1; direction  = 1;
function [value,isterminal,direction] = event1m(t,y,u,param)
  value      = y(1) - param.mulcopt;
  isterminal = 1; direction  = 1;
function [value,isterminal,direction] = event1r(t,y,u,param)
  value      = y(1) - param.coptrmd;
  isterminal = 1;direction  = 1;
function [value,isterminal,direction] = event1t(t,y,u,param)
  value      = y(1) - param.copt;
  isterminal = 1; direction  = 1;
function [value,isterminal,direction] = event2amt(t,y,u,param)
  value      = (y(1)/y(2)) - (param.c1f/param.c2f);
  isterminal = 1;
  direction  = 0;
function [value,isterminal,direction] = event2d(t,y,u,param)
  value      = (y(1)/y(2)) - (param.c1f/param.c2f);
  isterminal = 1; direction  = 0;
function [value,isterminal,direction] = event2m(t,y,u,param)
  value      = (y(1)/y(2)) - (param.c1f/param.c2f);
  isterminal = 1; direction  = 0;
function [value,isterminal,direction] = event2r(t,y,u,param)
  value      = (y(1)/y(2)) - (param.c1f/param.c2f);
  isterminal = 1; direction  = 0;
function [value,isterminal,direction] = event2t(t,y,u,param)
  value      = y(2) - param.c2f;
  isterminal = 1; direction  = 0;
function [value,isterminal,direction] = event3amt(t,y,u,param)
  value      = y(1) - param.c1f;
  isterminal = 1; direction  = -1;
function [value,isterminal,direction] = event3d(t,y,u,param)
  value      = y(1) - param.c1f;
  isterminal = 1; direction  = -1;
function [value,isterminal,direction] = event3m(t,y,u,param)
  value      = y(1) - param.c1f;
  isterminal = 1; direction  = -1;
function [value,isterminal,direction] = event3r(t,y,u,param)
  value      = y(1) - param.c1f;
  isterminal = 1; direction  = -1;
function [value,isterminal,direction] = event3t(t,y,u,param)
  value      = y(1) - param.c1f;
  isterminal = 1; direction  = 1;