Промышленный лизинг
Методички
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 |