% 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.li = ones(n,1)*(1/n)*6; 
optimparam.tf = 6;
optimparam.ui = 0.005*ones(1,n); 

optimparam.par = []; 
optimparam.bdu = [0 .5];
optimparam.bdx = [0 1600;0 600; 0.01 0.035];
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);

s= [68.1250e-9,   -56.4512e-6,   32.5553e-3,   -4.3529e-9, ...
    3.3216e-6,   -2.7141e-3];
s(6) = 1 * s(6);

q =  (s(1)*xplot(:,2).^2 + s(2)*xplot(:,2) + s(3)) .* exp((s(4)* ...
						  xplot(:,2).^2 + ...
						  s(5)*xplot(:,2) + ...
						  s(6)).*xplot(:,1)); 

alpha = uplot./q;

%..........................................................................
plot(tplot,alpha,'k')
title('')
xlabel('time')
ylabel('\alpha')
set(gcf,'PaperUnits','centimeters')
set(gcf,'PaperPosition',[0 0 11 8.25])
print -depsc2 a-alpha.eps
figure

plot(tplot,uplot,'k')
title('')
xlabel('time')
ylabel('u')
axis tight 
set(gcf,'PaperUnits','centimeters')
set(gcf,'PaperPosition',[0 0 11 8.25])
print -depsc2 a-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(gcf,'PaperUnits','centimeters')
set(gcf,'PaperPosition',[0 0 11 8.25])
print -depsc2 a-x12.eps
figure

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


