% Solve steady-state optimization problem for different disturbances to
% determine the active constraint regions as a function of disturbances
% Simulation results published in in Krishnamoorthy & Skogestad, Ind. Eng.
% Chem. Res., Bonvin Festschrift
%
% Author: Dinesh Krishnamoorthy, dineshk@ntnu.no
% Sep 2018, NTNU

clear
clc

addpath ('C:\Users\dineshk\CasADi\casadi-matlabR2014b-v3.1.0-rc1')
import casadi.*
d = 0.7:0.01:1.1;
for i = 1:numel(d)
    par.A1 = 5000;
par.A2 = 1e6;
cp = 1000;
par.E1 = 10000;
par.E2 = 15000;
R = 1.987;
Hr = 5000;
rho = 1;
V = 1*60;
pCB = 2.009;
pTi = 1.657e-3;
F_sp = 1;

par.Ts = 1;

dist = [d(i);0.0];

lbx = [0;0;0];
ubx = [100;100;1000];

lbu = [0;0;dist];
ubu = [1000;1;dist];

dx0 = [0.4978;0.5022;426.71];
u0 = [424.2;1;dist];

% x
CA = MX.sym('CA',1);
CB = MX.sym('CB',1);
T = MX.sym('T',1);

% dist
CAin = MX.sym('CAin',1);
CBin = MX.sym('CBin',1);

% u
Tin = MX.sym('Tin',1);
F = MX.sym('F',1);

% p
A1 = MX.sym('A1',1);
A2 = MX.sym('A2',1);

E1 = MX.sym('E1',1);
E2 = MX.sym('E2',1);

du = 1;
% ODE system
dx1 = F*(CAin - CA)/V - (A1*exp(-E1/(R*T))*CA - A2*exp(-E2/(R*T))*CB);
dx2 = F*(CBin - CB)/V + (A1*exp(-E1/(R*T))*CA - A2*exp(-E2/(R*T))*CB);
dx3 = F*(Tin - T)/V + Hr*((A1*exp(-E1/(R*T)))*CA - (A2*exp(-E2/(R*T)))*CB)/(rho*cp);

diff = vertcat(dx1,dx2,dx3);

diff = substitute(diff,A1,par.A1);
diff = substitute(diff,A2,par.A2);
diff = substitute(diff,E1,par.E1);
diff = substitute(diff,E2,par.E2);

x_var = vertcat(CA,CB,T);
p_var = vertcat(Tin,F,CAin,CBin);
d_var = vertcat(CAin,CBin);
L = -(pCB*CB - (pTi*Tin).^2) + 10000*(F - F_sp)^2;

ode = struct('x',x_var,'p',p_var,'ode',diff,'quad',L);  % der(m_tot) = w_in - w_out;
opts = struct('tf',par.Ts);

f = Function('f',{x_var,p_var},{diff,L},{'x','p'},{'xdot','qj'});

% create IDAS integrator
F_integrator = integrator('F','cvodes',ode,opts);


    
    %% Steady state optimization
    w = {};
    w0 = [];
    lbw = [];
    ubw = [];
    
    g = {};
    lbg = [];
    ubg = [];
    
    w = {w{:},x_var,p_var};
    lbw = [lbw;lbx;lbu];
    ubw = [ubw;ubx;ubu];
    w0 = [w0;dx0;u0];
    
    J = L; % Economic objective
    
    g = {g{:},vertcat(diff)};
    lbg = [0.*ones(3,1)];
    ubg = [0.*ones(3,1)];
    
    g = {g{:},x_var(3)};
    lbg = [lbg;0];
    ubg = [ubg;425];
    
    g = {g{:},x_var(1)};
    lbg = [lbg;0];
    ubg = [ubg;0.5];
    
    opts = struct('warn_initial_bounds',false, ...
        'print_time',false, ...
        'ipopt',struct('print_level',5) ...
        );
    
    nlp = struct('x',vertcat(w{:}),'f',J,'g',vertcat(g{:}));
    solver = nlpsol('solver','ipopt',nlp,opts);
    sol = solver('x0',w0,'lbx',lbw,'ubx',ubw,'lbg',lbg,'ubg',ubg);
    w_opt_SS = full(sol.x);
    J_opt = -full(sol.f);
    disp('-------Steady state Optimum-------')
    disp(' ')
    display([w_opt_SS(1:5);-full(sol.f)]);
    disp('----------------------------------')
    
    opt.w_opt(:,i) = w_opt_SS;
    opt.J(i) = -full(sol.f);
   end
opt.d = d;

save('CSTR_ACRegions.mat','opt')
%%

% Definition of figure and values
fig.dimw = 10;
fig.dimh = 6;
fig.Font = [16 14];
fig.label = [];
itg.fig = 1;
itg.lab = 1;
itg.arr = 1;
itg.fda = 1;

fig.in(itg.fig) = figure(itg.fig);
set(fig.in(itg.fig),'units','centimeters');
set(gcf, 'position', [5+itg.fig,5+itg.fig, fig.dimw, fig.dimh]);
clf
hold all
plot(opt.d,opt.w_opt(1,:)/0.5,'linewidth',1.5)
plot(opt.d,opt.w_opt(3,:)/425,'linewidth',1.5)
plot(opt.d,opt.w_opt(5,:),'k--','linewidth',1)
plot([0.85,0.85],[0,1.1],'k:')
plot([1,1],[0,1.1],'k:')
% plot(opt.d,opt.w_opt(1,:)-0.5,'linewidth',1.5)
% plot(opt.d,opt.w_opt(3,:)-425,'linewidth',1.5)
% plot(opt.d,opt.w_opt(5,:)-1,'k--','linewidth',1)
axs = gca;
axs.TickLabelInterpreter = 'latex';
% xticks([7:0.5:9])
ylabel('Normalized constraints','Interpreter','Latex')
xlabel('$C_{A,i}$','Interpreter','Latex')
% ylim([0.95,1.01])
ylim([0.65,1.05])
xlim([0.7,1.1])
leg = legend('$C_A/C_A^{max}$','$T/T^{max}$','$F/F^{sp}$');
set(leg,'Interpreter','latex','Orientation','Vertical','Location','Best','box','off')
box on

