close all
clc
clear all
%Try to optimize distillation column "A" for different disturbances: Here
%we use a different objective function.

% Define nominal inputs!

% Nominal inputs
LT=2.70629;                          % Reflux
VB=3.20629;                          % Boilup
D=0.5;                               % Distillate
B=0.5;                               % Bottoms
F=1.0;                               % Feedrate
zF=0.5;                              % Feed composition
qF=1.0;                              % Feed liquid fraction
Uinit = [LT VB D B F zF qF]';
global NT d
NT = 41 ;

%Now we split this set in three: LT and VB are our unused inputs which can
%be used for optimization. D and B are used to control the holdups in the
%column. F, zF and qF are considered disturbances (but we leave qF
%unchanged for now).

u0 = Uinit(1:2) ;

options = optimset('Display','off','TolFun',1e-6,'TolCon',1e-6,'Algorithm','active-set','MaxIter',75) ;

lb = [0.1 0.1]' ;
ub = [10 u0(2)*1.25]' ;

pV = 0.12;
F = 1.0 ;

act = zeros(length(pV),length(F)) ;
kmax = length(pV)*length(F) ;
k = 0 ;

if length(pV) == 1 && length(F) > 1
    uopt = zeros(2,length(F)) ;
    copt = zeros(2,length(F)) ;
    lopt = zeros(3,length(F)) ;
elseif length(F) == 1 && length(pV) > 1
    uopt = zeros(2,length(pV)) ;
    copt = zeros(2,length(pV)) ;
    lopt = zeros(3,length(pV)) ;
end

if length(pV) == 1 && length(F) == 1
    options.Display = 'iter' ;
end

%Here the difference shows. We need to make room for more regions, as there
%is now a possibility that xd>xd_min is also inactive

for i=1:length(pV)
    for j=1:length(F)
    d = [F(j) zF qF pV(i)]' ;
    k = k+1;
    [u_opt,FVAL,EXITFLAG,OUTPUT,LAMBDA] = fmincon(@col_obj_II,u0,[],[],[],[],lb,ub,@col_constr,options);
        if EXITFLAG <= 0 
        disp(['Optimization number ' num2str(k) ' did not converge in ' num2str(options.MaxIter) ' iterations!']) 
        act(i,j) = NaN ;
        elseif EXITFLAG > 0
        disp(['Optimization number ' num2str(k) ' of ' num2str(kmax) ' completed']) 
        %Now check which constraint is active
            if LAMBDA.ineqnonlin(1) == 0 && LAMBDA.ineqnonlin(2) == 0 && LAMBDA.upper(2) == 0 % No active constraints 
                act(i,j) = 1 ;
            elseif LAMBDA.ineqnonlin(1) == 0 && LAMBDA.ineqnonlin(2) == 0 && LAMBDA.upper(2) > 0 % Vmax active
                act(i,j) = 2 ;
            elseif LAMBDA.ineqnonlin(1) == 0 && LAMBDA.ineqnonlin(2) > 0 && LAMBDA.upper(2) > 0 % Vmax and XD active
                act(i,j) = 3 ;
            elseif LAMBDA.ineqnonlin(1) == 0 && LAMBDA.ineqnonlin(2) > 0 && LAMBDA.upper(2) == 0 % XD active
                act(i,j) = 4 ;
            elseif LAMBDA.ineqnonlin(1) > 0 && LAMBDA.ineqnonlin(2) > 0 && LAMBDA.upper(2) == 0 % XD and XB active
                act(i,j) = 5 ;
            elseif LAMBDA.ineqnonlin(1) > 0 && LAMBDA.ineqnonlin(2) > 0 && LAMBDA.upper(2) > 0 % All constraints active! 
                act(i,j) = 6 ;
            elseif LAMBDA.ineqnonlin(1) > 0 && LAMBDA.ineqnonlin(2) == 0 && LAMBDA.upper(2) > 0 % Vmax and XB active 
                act(i,j) = 7 ;
            elseif LAMBDA.ineqnonlin(1) > 0 && LAMBDA.ineqnonlin(2) == 0 && LAMBDA.upper(2) == 0 % XB active
                act(i,j) = 8 ; 
            end    
        end  
        
        if length(pV) == 1 && length(F) > 1
                uopt(:,j) = u_opt ;
                C = col_constr(u_opt) ;
                copt(:,j) = C ;
                lopt(:,j) = [LAMBDA.upper(2) ; LAMBDA.ineqnonlin] ;
        elseif length(F) == 1 && length(pV) > 1
                uopt(:,i) = u_opt ;
                C = col_constr(u_opt) ;
                copt(:,i) = C ;
                lopt(:,i) = [LAMBDA.upper(2) ; LAMBDA.ineqnonlin] ;
         end
     end
end

if length(pV) >1 && length(F) > 1
act = flipud(act) ;
disp(act)
elseif length(pV) == 1    && length(F) > 1
figure(1),plot(F,copt(1,:),'r',F,copt(2,:),'b') ;
figure(2),plot(F,uopt(1,:),'r',F,uopt(2,:),'b') ;
figure(3),plot(F,lopt(1,:),'r',F,lopt(2,:),'b',F,lopt(3,:),'gr') ;
matris = [act ; lopt]';
disp(matris)
elseif length(F) == 1 && length(pV) > 1
figure(1),plot(pV,copt(1,:),'r',pV,copt(2,:),'b') ;
figure(2),plot(pV,uopt(1,:),'r',pV,uopt(2,:),'b') ;  
figure(3),plot(pV,lopt(1,:),'r',pV,lopt(2,:),'b',pV,lopt(3,:),'gr') ;
act = act' ;
matris = [act ; lopt] ;
disp(matris)

elseif length(F) == 1 && length(pV) == 1
[Jcol D B xD xB x] = col_obj_II(u_opt) ;
end


