% Needs dynopt optimisation package
options = optimset('LargeScale','off','Display','iter');
options = optimset(options,'GradObj','on','GradConstr','on');
options = optimset(options,'MaxFunEvals',1e5);
options = optimset(options,'MaxIter',1e5);
options = optimset(options,'TolFun',1e-7);
options = optimset(options,'TolCon',1e-7);
options = optimset(options,'TolX',1e-7);
options = optimset(options,'Algorithm','active-set'); %2008b

n=3;
optimparam.optvar = 3; 
optimparam.objtype = []; 
optimparam.ncolx = 3; 
optimparam.ncolu = 2; 
optimparam.tf = [];
%optimparam.li = ones(n,1)*(1/n)*4.5; 
%optimparam.ui = 0.005*ones(1,n); 
optimparam.li =[1.6 1.4 1.5 ]';
optimparam.ui=[0 0.009 0.009];

optimparam.par = []; 
optimparam.bdu = [0 .5];
optimparam.bdx = [0 1600;0 600; 0.01 0.035; 0 10];
optimparam.bdp =[];
optimparam.objfun = @objfun;
optimparam.confun = @confun;
optimparam.process = @process;
optimparam.options = options;


[optimout,optimparam]=dynopt(optimparam)
[tplot,uplot,xplot] = profiles(optimout,optimparam,50);
%[tp,cp,ceqp] = constraints(optimout,optimparam,50);

%..........................................................................
plot(tplot,uplot,'k')
title('')
xlabel('time')
ylabel('u')
axis tight 
set(1,'PaperUnits','centimeters')
set(1,'PaperPosition',[3 10 11 8.25])
print -depsc2 b-control.eps

figure
plot(tplot,xplot(:,1),'b:',tplot,xplot(:,2),'r--')
title('')
xlabel('time')
ylabel('x_1, x_2')
legend('x_1','x_2',4)
axis tight 
set(1,'PaperUnits','centimeters')
set(1,'PaperPosition',[3 10 11 8.25])
print -depsc2 b-x12.eps

figure
plot(tplot,xplot(:,3),'k-')
title('')
xlabel('time')
ylabel('x_3')
set(1,'PaperUnits','centimeters')
set(1,'PaperPosition',[3 10 11 8.25])
v=axis;
axis([v(1) v(2) 0.0098 v(4)])
print -depsc2 b-x3.eps

%..........................................................................
