Промышленный лизинг Промышленный лизинг  Методички 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 [ 62 ] 63 64 65 66 67

plot(State(:,l),State(:,2),x-) xlabel(State l) ylabel (State 2 )

figure

t2=[0:0.001:0.999];

plot(t2,Control(floor(nn*t2)+1,1),r) hold on

plot(t2,Control(floor(nn*t2)+1,2),g) %The following codes comprise the SCOM package; the user does not alter them

% Get function values as required by constr

function [f,g]=fqq(uz,ps,q,xinit) % Input data to the subroutine if ps(5)==1

[a1,a2,a3,a4,a5,a6]=cqq(uz,ps,q,xinit); % Output data elseif ps(5)==0

[a1,a2,a3,a4]=cqq(uz,ps,q,xinit); end

f=a1;

g=a2;

function [df,dg]=gqq(uu,ps,q,xinit)

% Get gradient values as required by constr

[a1,a2,a3,a4,a5,a6]=cqq(uu,ps,q,xinit); df=a6; nn=ps(3); pk=char(q(5)); dg=feval(pk,uu,nn);

function [f,g,xm,f0,lm,gr]=cqq(um,ps,q,xinit) % Solve differential equations

nx=ps(1);nu=ps(2);nn=ps(3);npa=ps(4); rec=0;npar=1; xm=zeros(nn+1,nx); xm(1,:)=xinit; lm=zeros(nn+1,nx); ma=nx;t=0;it=1;hs=1/nn; px=char(q(1)) ;

xm=nqq(px,nx,nu,nn,xm,ma,t,it,hs,um,xm,lm,ps); %compute state

ma=1;t=0;it=1;hs=1/nn;

zz=zeros(nn+1,1); zz(1,:)=0;



pj=char (q(2));

jm=nqq(pj,nx,nu,nn,zz,ma,t,it,hs,um,xm,lm,ps);%compute integral

xf=xm(nn+1,:) ;

pf=char(q(3));

f=jm(nn+1) + feval(pf,xf,um,xm,ps); %objective pc=char(q(4));

for ii=1:nn

g(ii)=feval(pc,ii,hs,um,xm,lm,ps) ; %control constraint end

f0=jm(nn+1); %final state

if ps(5)==1

%if gradients are supplied

ma=nx;t=1;it=nn;hs=-1/nn;

lm=zeros(nn+1,nx);

pa=char(q(8));

lm(nn+1,:)=feval(pa,nn,xf,um,xm,ps); pq=char(q(6));

lm=lqq(pq,nx,nu,nn,lm,ma,t,it,hs,um,xm,lm,ps); %costate

function xk=nqq(pd,nx,nu,nn,fil,ma,t,it,hs,um,xm,lm,ps) %Runge-Kutta (organize steps for fourth order RK integration of DEs)

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 : increments

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) % Runge-Kutta : get function values

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 xk=lqq(pd,nx,nu,nn,fil,ma,t,it,hs,um,xm,lm,ps) %Runge-Kutta for adjoint differential equation (time reversed)

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

while it > 0 %< nn+1

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

function ix=iqq(vv,hh,xi)

% Linear interpolation (forwards)

nn=length(vv)-1;

nz=1/hh;

fr=nz*mod(xi,hh); bb=floor(xi*nz)+1;

if bb <= nn, vu=vv(bb+1,:);

else vu=vv(nn);

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

function ix=jqq(vv,hh,xi)

% Linear interpolation (backwards)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 [ 62 ] 63 64 65 66 67