動(dòng)態(tài)規(guī)劃matlab仿真實(shí)例

1、動(dòng)態(tài)規(guī)劃在火力分配中的應(yīng)用。1.問題描述設(shè)有m個(gè)目標(biāo)，目標(biāo)價(jià)值(重要性和危害性)各不相同，用數(shù)值A(chǔ)k(K=1,2,.m)表示，計(jì)劃用n枚導(dǎo)彈突襲，導(dǎo)彈擊毀目標(biāo)的概率PK=1-e-%叫，其中玉是常數(shù)，取決于導(dǎo)彈的特性與目標(biāo)的性質(zhì)；叫為向目標(biāo)發(fā)射的導(dǎo)彈數(shù)，問題：做出方案使預(yù)期的突擊效果最大。2. 問題建模上述問題可以表述為maxV=Y化'JI一e-叫叫)約束條件為=n(叫為非負(fù)整數(shù))3. 算法描述下面通過一個(gè)實(shí)例說明：設(shè)目標(biāo)數(shù)目為4(m=4)，導(dǎo)彈為5(n=5)，4上和aK取值情況如下表所示：表1：Ak，叫取值情況目標(biāo)k123487630.20.30.50.9將火力分配可分為4個(gè)階段，每個(gè)

2、階段指標(biāo)函數(shù)為：片("J=8(1一£-0.2“1)2(U2)=7(1£-0.3“2)匕(出)=6(1一£-0.5“3)/(氣)=3(1一60-94)叫可能取值為0，1,2,3,4,5,將函數(shù)值帶人如下表:表2函數(shù)值u叫(叫)晌2)匕4)伽)0000011.451.812.361.7922.643.163.792.5133.614.154.662.8144.414.895.192.9355.065.445.512.97動(dòng)態(tài)規(guī)劃問題基本方程為：f(呂)=maxV(u+f(x-u)ckkkkk+1kkf(%J=055逐次向前推一級(jí)K=4小4)=V4(u4)=3

3、(1-曠0叫)K=3f(%)=maxV(u)+f(x333343)=max6(1e0-5u3)+/(xu)433K=2f(x)=maxV(u)+f(x222232)=max7(1e0-3u2)+/(xu)322K=1f(x)=maxV(u>)+f(x-111121-%)=max8(1e0-2u1)+/(xu)211(0<u<5可取等號(hào))k只需要求解人(5)的最大值然后反推回去就可以獲得最優(yōu)的分配方案4.Matlab仿真求解因?yàn)榻信c叫取值為整數(shù)，可以采用動(dòng)態(tài)規(guī)劃的方法，獲得人(5)的最大值，對(duì)應(yīng)的最優(yōu)方案functionp_opt,fval=dynprog(x,DecisFun

4、,SubObjFun,TransFun,ObjFun)%求解動(dòng)態(tài)規(guī)劃問題最小值函數(shù)k=length(x(1,:)%判斷決策級(jí)數(shù)x_isnan=isnan(x);%非空狀態(tài)矩陣t_vubm=inf*ones(size(x);%性能指標(biāo)中間矩陣f_opt=nan*ones(size(x);%總性能指標(biāo)矩陣d_opt=f_opt;%每步?jīng)Q策矩陣tmp1=find(x_isnan(:,k);%最后一步狀態(tài)向量tmp2=length(tmp1);%最后一步狀態(tài)個(gè)數(shù)fori=1:tmp2u=feval(DecisFun,k,x(tmp1(i),k);tmp3=length(u);決策變量forj=1:tm

5、p3%求出當(dāng)前狀態(tài)下所有決策的最小性能指標(biāo)tmp=feval(SubObjFun,k,x(tmp1(i),k),u(j);iftmp<=t_vubm(i,k)%t_vubf_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vubm(i,k)=tmp;end;end;endforii=k-1:-1:1tmp10=find(x_isnan(:,ii);tmp20=length(tmp10);fori=1:tmp20%求出當(dāng)前狀態(tài)下所有可能的決策u=feval(DecisFun,ii,x(tmp10(i),ii);tmp30=length(u);forj=1:tmp30%求出當(dāng)前

6、狀態(tài)下所有決策的最小性能指標(biāo)tmp00=feval(SubObjFun,ii,x(tmp10(i),ii),u(j);%單步性能指標(biāo)tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j);下%一狀態(tài)tmp50=x(:,ii+1)-tmp40;找%出下一狀態(tài)在x矩陣的位置tmp60=find(tmp50=0);ifisempty(tmp60)ifnargin<6矩陣不同需要修改nargin的值，很重要tmp00=tmp00+f_opt(tmp60(1),ii+1);%setthedefaultobjectvalueelsetmp00=feval(ObjFun

7、,tmp00,f_opt(tmp60(1),ii+1);end當(dāng)前%狀態(tài)的性能指標(biāo)iftmp00<=t_vubm(i,ii)f_opt(i,ii)=tmp00;d_opt(i,ii)=u(j);t_vubm(i,ii)=tmp00;end;end;end;end;endfval=f_opt(:,1);tmp0=find(isnan(fval);fval=fval(tmp0,1);p_opt=;tmpx=;tmpd=;tmpf=;tmp01=length(tmp0);fori=1:tmp01tmpd(i)=d_opt(tmp0(i),1);tmpx(i)=x(tmp0(i),1);tmpf

8、(i)=feval(SubObjFun,1,tmpx(i),tmpd(i);p_opt(k*(i-1)+1,1,2,3,4)=1,tmpx(i),tmpd(i),tmpf(i);forii=2:ktmpx(i)=feval(TransFun,ii,tmpx(i),tmpd(i);tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1=0);ifisempty(tmp2)tmpd(i)=d_opt(tmp2(1),ii);endtmpf(i)=feval(SubObjFun,ii,tmpx(i),tmpd(i);p_opt(k*(i-1)+ii,1,2,3,4)=ii,tmpx(

9、i),tmpd(i),tmpf(i);end;end;下面編寫四個(gè)函數(shù)functionu=DecisF1(k,x)%決策函數(shù)ifk=4u=x;elseu=0:x;%狀態(tài)轉(zhuǎn)移方程%階段k的指標(biāo)函數(shù)基%本方程中的函數(shù)endfunctiony=TransF1(k,x,u)y=x-u;functionv=SubObjF1(k,x,u)a=0.2,0.3,0.5,0.9;A=8,7,6,3;v=A(k)*(1-exp(-a(k)*u);v=-v;%max變?yōu)閙infunctiony=ObjF1(v,f)y=v+f;y=-y;%max變?yōu)閙in測(cè)試代碼：clear;n=5;x1=n;nan*ones(n,1);x2=0:n;x2=x2'x=x1,x2,x2,x2;p,f=dynprog(x,'DecisF1','SubObjF1','TransF1','ObjF1')%p為分配方案，f為結(jié)果5.Matlab仿真結(jié)果分析運(yùn)行結(jié)果顯示：P為方案：田pu4k4doubS1|1234151-145022241-1.81433332-3.79274411-1.7803即對(duì)目標(biāo)1發(fā)射1枚導(dǎo)彈，對(duì)目標(biāo)2發(fā)射1枚，對(duì)目標(biāo)3發(fā)射2枚，對(duì)目標(biāo)4發(fā)生1枚導(dǎo)彈，能獲得