% feas_script.m
% This is version v1.0
% You need Tomlab installed in order to run this script

% Changelog:
global R d vub vlb hessflag teller coldata d cs_glob unc_glob R_glob


% Load data in memory
startup % Define the parameters to use in this script
disp('Calculating the nominal optimum')
% GLOBAL FLAG
flagrecalc=0; % 0-Recalculate everything, 1- load from file
% Find the nominal optimum
[znomopt,Vnomopt,lambda,XA]=nom_opt(1); % flag=1, print output
xnomopt=znomopt(1:104);
unomopt=znomopt(105:113);
%save xopt xnomopt
%save uopt unomopt
%else
%    load xopt
%    load uopt
%    znomopt=[xnomopt;unomopt;coldata.yset];
%end
% START LOSS CALCULATION
% Vi kjører test på to strukturer R(1,2) og R(3,4) hvor den første er DTS
% og den andre er å holde Rl og Rv konstant
R=zeros(10,60);
indDTS=[4 12 28 36];
Temp(4)+Temp(12)-Temp(28)-Temp(36);
indDTS2=[13 28 32 41];
nIND=[1;1;-1;-1];
nIND2=[1;-1;1;-1];
R(1,3+indDTS)=nIND';         % DTS
R(2,56)=1;                   % Rv
R(3,53)=1;                   % V_53
R(4,55)=1;                   % Rl_55
R(5,3+indDTS2)=nIND2';

% IMPLEMENTATION ERROR DIFFERENT STRUCTURES
% Disse er ikke aktive (ganger med null foran)
imperr=zeros(4,3);
imperr(1:2,:)=0*[-0.2 0 0.2;-0.2 0 0.2];
imperr(3:4,:)=0*[-0.2 0 0.2;-0.02 0 0.02];
disp('Calculate loss')
%%%%%%% Disturbances
% Her definerer man forstyrrelser.Her har jeg antatt at man har 2
% forstyrrelser som kan skje simultant.I dette tilfelle 
% dvec1=[0;0;-0.1;0;0;0;0]; % Forstyrrlese i ql
%      F zA zB ql xAset xBset xCset [+0.1;+0.1;-0.1;-0.1;0;0;0]

% Nominal measurements
Ynom=genpet(0,xnomopt,unomopt,3,coldata.Nt,coldata.alfa,1,1,[]); % Y=[xA xB xC Temp u]=[1 1 1 48 9]
samle=[];
for i=0:0
hessflag=i;
% Optimization parameters
if i==0
c_L = zeros(114,1); % Lower bound inequality constraints
c_U = zeros(114,1);      % Upper bound inequality constraints
iii=1;
disp('======================================================================')
disp('                  Ingen beskrankning på pådrag                        ')
disp('======================================================================')
elseif i==1
c_L = zeros(116,1); % Lower bound inequality constraints
c_U = zeros(116,1);      % Upper bound inequality constraints
iii=1;
disp('======================================================================')
disp('                        Fixed Rv + DTS                                ')
disp('======================================================================')
elseif i==2
c_L = zeros(116,1); % Lower bound inequality constraints
c_U = zeros(116,1);      % Upper bound inequality constraints
iii=1;
disp('======================================================================')
disp('                        Fixed  Rv + Rl                                 ')
disp('======================================================================')
elseif i==3
c_L = zeros(116,1); % Lower bound inequality constraints
c_U = zeros(116,1);      % Upper bound inequality constraints
iii=1;
disp('======================================================================')
disp('                        Fixed  Rv + DTS2                              ')
disp('======================================================================')
elseif i==4
c_L = zeros(116,1); % Lower bound inequality constraints
c_U = zeros(116,1);      % Upper bound inequality constraints
iii=20;
disp('======================================================================')
disp('                        Fixed  Rv + V                                 ')
disp('======================================================================')
end

% c_U(112:114)=0.03;
lb=zeros(116,1);
ub=10*ones(116,1);
Name = 'Loss calculations';
%Prob = conAssign('petl_opt',[],[],[],lb,ub,Name,zopt,[],[],[],[],[],'petl_cons',[],[],[],c_L,c_U,[],[],[],[]);
Prob.PriLevOpt = 2;
Prob.PriLev=2;
Prob.Solver.Alg =4;
Prob.NumDiff =6; % Use standard numerical differences
Prob.ConsDiff=6; % Use complex variable method to estimate differences
Prob.SOL.optpar = -999*ones(62,1);
Prob.SOL.optpar(1)=111111;
Prob.SOL.PrintFile ='opt.pri';
Prob.SOL.optPar(39)=0;
Prob.SOL.optPar([9,11,12])=1e-6;
eflagc=1;ii=1;dv=1;
            while eflagc>0 & ii<=iii
%                 dv=1+(ii-1)*0.05;
            cs_glob(1)=R(1,:)*Ynom;      % DTS
            cs_glob(2)=R(2,:)*Ynom;      % Rv
            cs_glob(3)=R(3,:)*Ynom;      % V
            cs_glob(4)=R(4,:)*Ynom;      % Rl
            cs_glob(5)=R(5,:)*Ynom;      %DTS2
            R_glob=R;
            unc_glob(1)=0;%usetvec1(l3);
            unc_glob(2)=0;%usetvec2(l4);
            d=[unomopt(6:9);XA;coldata.yset(2:3)]+dvec1(:,1);
            % fmincon
            
            Prob = conAssign('petl_opt',[],[],[],lb,ub,Name,znomopt,[],[],[],[],[],'petl_cons',[],[],[],c_L,c_U,[],[],[],[]);
            Result_Loss=tomRun('snopt',Prob);
%             textfil=num2str(globalteller);
%            filnavn=['Resultat' textfil]
%            cd ./Results
%            save(filnavn,'Result_Loss')
%            cd ..
            zoptc(:,1)=Result_Loss.x_k;
            uoptc(:,1)=Result_Loss.x_k(105:113);
            optc=Result_Loss.f_k;
            eflagc=Result_Loss.ExitFlag;
            Result_Loss.ExitFlag;
            disp(Result_Loss.ExitText);
            clear Result_Loss
            ii=ii+1;
        end
disp('eflagc(1)=0 OK, if not not converged or infeasible, check Result_Loss.ExitText')

eflagc

uoptc=zoptc(105:113);
L0=uoptc(1);			% Reflunomoptx rate 
VB0=uoptc(2);			% Boilup rate 
S0=uoptc(3);			% Sideproduct removal rate 
Rl0=uoptc(4);			% Liquid split fraction to feed tray side 
Rv0=uoptc(5);  		% Vapour split fraction to feed tray side 
F0=uoptc(6);			% Feed flow rate
zF0=uoptc(7:8);	    % Feed composition for components 1:nc-1
qF0=uoptc(9);	
disp(sprintf('L %s V %s S %s Rl %s Rv %s F %s zF1 %s zF2 %s qF %s dv',9,9,9,9,9,9,9,9,9));
disp(sprintf('%6.3g %6.3g %6.3g %6.3g %6.3g %6.3g %s %6.3g %6.3g %6.3g',L0,VB0,S0,Rl0,Rv0,F0,9,zF0(1),zF0(2),qF0,dv));
disp(sprintf('xA %s xC %s xB',9,9));
disp(sprintf('%6.3g %6.3g %6.3g',zoptc(114),zoptc(115),zoptc(116)));
samle=[samle;L0,VB0,S0,Rl0,Rv0,F0,zF0(1),zF0(2),qF0,zoptc(114),zoptc(115),zoptc(116),eflagc,dv];
end