Appendix C

Differential Equation Solver from the SCOM Package

function xk=nqq(pd,nx,nu,nn,fil,ma,t,it,hs,um,xm,lm,ps) % Calling the Runge-Kutta integration

yin=fil(1,:);xk(1,:)=yin; while it < nn+1

[y2,it2,t2]=rqq(pd,ma,t,it,hs,yin,nx,nu,nn,um,xm,lm,ps); xk(it+1,:)=y2;

it=it2;t=t2;

yin=y2;

function [y1,it,t]=rqq(pd,ma,t,it,hs,yin,nx,nu,nn,um,xm,lm,ps) % Runge-Kutta stages

fp=zeros(ma,1); tt=zeros(ma,1); P=0; q=1;

tz=kqq(pd,ma,t,it,hs,p,q,fp,yin,tt,um,xm,lm,ps);

tt=tz(1,:); fp=tz(2,:);

p=0.5;q=2;t=t+0.5*hs;

tz=kqq(pd,ma,t,it,hs,p,q,fp,yin,tt,um,xm,lm,ps); tt=tz(1,:); fp=tz(2,:);

tz=kqq(pd,ma,t,it,hs,p,q,fp,yin,tt,um,xm,lm,ps);

tt=tz(1,:); fp=tz(2,:);

t=t+0.5*hs;p=1;q=1;

tz=kqq(pd,ma,t,it,hs,p,q,fp,yin,tt,um,xm,lm,ps);

tt=tz(1,:); fp=tz(2,:);

it=it+sign(hs);

y1=yin+tt/6;

function tz=kqq(pd,ma,t,it,hs,p,q,fp,yin,tt,um,xm,lm,ps) % Calling the function

z=yin+p*hs*fp;

ff=feval(pd,t,it,z,yin,hs,um,xm,lm,ps); fp=ff;

tz=[tt+q*hs*fp;fp] ;

function ix=jqq(vv,hh,xi) % Linear interpolation

global um xm jm

nn=length(vv)-1;

nz=-1/hh;

if xi < 0, xi = 0; end

fr=-nz*mod(xi,hh);

bb=ceil(xi*nz)+1;

if bb>1, vw=vv(bb-1,:);

else vw=vv(1);

ix=(1-fr)*vv(bb,:) + fr*vw;

Appendix D SCOM Package

subs=cell(1,9);

subs={dnx , dnj , dnf , dnc};

par=[2, 2, 20, 0, 0, 1, 6, 0.2, 0.1, 0.1, 0.75]; % nx, nu, nn, npa, grad, c, T, r, p, d , b=k/r\ % nx = number of states, nu = number of controls,

nn = number of subdivisions % npa = 0 (reserved), grad = 1 if gradients calculated, otherwise 0

% The remaining parameters, specific to the problem, are passed to the subroutines

xinit=[2 1]; nn=par(3); % Initial values for the state(s) u0=zeros(nn,2); % Starting values for computing the control ul=zeros(nn,2) ; % Lower bound (vector or matrix) for the control uu=ones(nn,2); % Upper bound (vector or matrix) for the control

figure

Control=constr(fqq ,u0, [] ,ul,uu, [] ,par,subs,xinit) % Calls constr package

[Objective,Constraint,State,Integral]= cqq(Control,par,subs,xinit) % Computes optimal state, etc., from the optimal control got from constr

% If gradients are calculated, then used instead

% [Objective, Constraint, State, Integral, Costate,Gradient]

=cqq(Control,par,subs,xinit) % The following lines plot one state component against another,

and plot two control components, ehich are step-functions,

against time