%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Steady-state optimization - region

clear all
global Y g

if 1==0
cd /home/msg/Matlab
startup
cd /home/msg/Matlab/larsson2002/caseC
pwd
end

%options
options=optimset;
options.DerivativeCheck='off';
options.Diagnostics='on';
options.Display='iter';
options.GradConstr='on';
options.GradObj='on';
options.LargeScale='off';
options.TolCon=1e-9;
options.TolFun=1e-9;
options.TolX=1e-9;

%Disturbances - regions
%p1=1400; p1=p1(:);
%p20=[490:1:495]; p20=p20(:);
%p1=[1399.5:0.1:1400.5]; p1=p1(:);
%p20=[492.6:0.01:492.7]; p20=p20(:);
%p1=[399.5:0.1:400.5]; p1=p1(:);
%p20=[242.1:0.01:242.2]; p20=p20(:);
%p1=[1399.5:0.1:1400.5]; p1=p1(:);
%p20=[487.9:0.01:488.0]; p20=p20(:);
%%p1=[399.5:0.1:400.5]; p1=p1(:);
%%p20=[242.1:0.01:242.2]; p20=p20(:);

p1=1400; p1=p1(:);
p20=[490:-1:460]; p20=p20(:);
%p20=[463:-0.1:462]; p20=p20(:);
%p20=[487:0.1:488]; p20=p20(:);

%Initial values
load results/opt0
F0max=492;
p=p0; p(20)=F0max;%Nominal
ux=ux0;
np1=size(p1,1);
np20=size(p20,1);
YY=[];

%Optimization
tcpu=cputime;
for i=1:np1
   p(1)=p1(i);
   for j=1:np20
      p(20)=p20(j);
      [ux,fval,exitflag,output,lambda,grad,hessian]=fmincon('fun',ux,[],[],[],[],[],[],'nonlcon',options,p);
      [ux,fval,exitflag,output,lambda,grad,hessian]=fmincon('fun',ux,[],[],[],[],[],[],'nonlcon',options,p);
      VBmax=p(1); F0max=p(20); Mr=Y(7); xb=Y(12); V=Y(2); F0=Y(6);
      disp(sprintf('VBmax=%4.2f F0max=%4.2f Mr=%6.4f xb=%6.4f V=%6.4f F0=%6.4f exitflag=%d max(g)=%6.4f',VBmax,F0max,Mr,xb,V,F0,exitflag,max(g(3:8))));   
      YY=[YY;VBmax F0max Mr xb V F0 exitflag max(g(3:8))]; 
   end
end
cputime-tcpu

%save results/opt_V1400_F492 p1 p20 YY 
%save results/opt_V400_F242 p1 p20 YY 

disp(sprintf('Vmax Fmax Mr xb V F0 exitflag max(g(3:8))'))
for i=1:size(YY,1)
   disp(sprintf('%9.4f %6.4f %6.4f %6.4f %9.4f %9.4f %6.4f %9.4f',YY(i,:)));
end  






