%% Main program
clear all; clc

% Inputs and disturbances
LT = 2.70629;       % Reflux
VB = 3.20629;       % Boilup
F  = 1.0;           % Feedrate
qF = 1.0;           % Feed liquid fraction
zF = [0 0.5 0.5 0]; % Feed compositions for component x1 to xNC
u1 = [LT VB];
d  = [F qF zF];

load steadystate % y0

%% Steady-state calculations
options = optimset('Display','iter','Diagnostics','on');
x  = fmincon(@cost,y0,[],[],[],[],[],[],@multicomA,options,u1,d);
f1 = output(x,u1,d);
disp('       L         V         D         B                 xD                            xB   ')
disp([f1{1}, f1{2}, f1{3}, f1{4}, f1{5}, f1{6}])

%% Optimization calculations
lb = zeros(length([x;u1']));
y  = fmincon(@cost,[x;u1'],[],[],[],[],lb,[],@multicomA,options,d);
z0 = y(1:length(y0)); u2 = y(length(y0)+1:end);
f2 = output(z0,u2,d);
disp('       L         V         D         B                 xD                            xB   ')
disp([f2{1}, f2{2}, f2{3}, f2{4}, f2{5}, f2{6}])

%% Dynamic calculations
[t,z] = ode15s(@multicomA,[0 100],x,[],u1,d);

