動(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-ak uk，其中ak是常數(shù)，取決于導(dǎo)彈的特性與目標(biāo)的性質(zhì)；uk為向目標(biāo)發(fā)射的導(dǎo)彈數(shù)，問題：做出方案使預(yù)期的突擊效果最大。2. 問題建模上述問題可以表述為 maxV=k=1mAk(1-e-ak uk)約束條件為 k=1muk=n （uk 為非負(fù)整數(shù)）3. 算法描述下面通過一個(gè)實(shí)例說明：設(shè)目標(biāo)數(shù)目為4（m=4），導(dǎo)彈為5（n=5），Ak和aK取值情況如下表所示： 表1：Ak ，uk取值情況目標(biāo)K1234Ak8763uk0

2、.20.30.50.9將火力分配可分為4個(gè)階段，每個(gè)階段指標(biāo)函數(shù)為：V1u1=8(1-e-0.2u1) V2u2=7(1-e-0.3u2) V3u3=6(1-e-0.5u3) V4u4=3(1-e-0.9u4)1 / 6uk 可能取值為0，1，2，3，4，5，將函數(shù)值帶人如下表：表2 函數(shù)值uV1u1V2u2V3u3V4u40000011.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ī)劃問題基本方程為：fkxk=maxVkuk+fk+1xk-uk cf5x5=0逐次向前推

3、一級(jí)K=4 f4x4=V4u4=3(1-e-0.9u4) K=3 f3x3=maxV3u3+f4x3-u3=max6(1-e-0.5u3)+ f4x3-u3K=2 f2x2=maxV2u2+f3x2-u2= max7(1-e-0.3u2)+ f3x2-u2K=1 f1x1= maxV1u1+f2x1-u1= max8(1-e-0.2u1)+ f2x1-u1 （0<uk<5可取等號(hào)）只需要求解f15的最大值然后反推回去就可以獲得最優(yōu)的分配方案4. Matlab仿真求解因?yàn)閤k 與uk取值為整數(shù)，可以采用動(dòng)態(tài)規(guī)劃的方法，獲得f15的最大值，對(duì)應(yīng)的最優(yōu)方案functionp_opt,fv

4、al=dynprog(x,DecisFun,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ù)for i=1:tmp2 u=feval(DecisFun,k,x(tmp

5、1(i),k); tmp3=length(u);%決策變量 for j=1:tmp3 % 求出當(dāng)前狀態(tài)下所有決策的最小性能指標(biāo) tmp=feval(SubObjFun,k,x(tmp1(i),k),u(j); if tmp <= t_vubm(i,k) %t_vub f_opt(i,k)=tmp; d_opt(i,k)=u(j); t_vubm(i,k)=tmp; end; end;end for ii=k-1:-1:1 tmp10=find(x_isnan(:,ii); tmp20=length(tmp10); for i=1:tmp20 %求出當(dāng)前狀態(tài)下所有可能的決策 u=feval

6、(DecisFun,ii,x(tmp10(i),ii); tmp30=length(u) ; for j=1:tmp30 % 求出當(dāng)前狀態(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) if nargin<6 %矩陣不同需要修改nargin的

7、值，很重要 tmp00=tmp00+f_opt(tmp60(1),ii+1); % set the default object value else tmp00=feval(ObjFun,tmp00,f_opt(tmp60(1),ii+1); end %當(dāng)前狀態(tài)的性能指標(biāo) if tmp00<=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)

8、;p_opt=;tmpx=;tmpd=;tmpf=; tmp01=length(tmp0); for i=1:tmp01 tmpd(i)=d_opt(tmp0(i),1); tmpx(i)=x(tmp0(i),1); tmpf(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); for ii=2:k tmpx(i)=feval(TransFun,ii,tmpx(i),tmpd(i); tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1=0); if

9、 isempty(tmp2) tmpd(i)=d_opt(tmp2(1),ii); end tmpf(i)=feval(SubObjFun,ii,tmpx(i),tmpd(i); p_opt(k*(i-1)+ii,1,2,3,4)=ii,tmpx(i),tmpd(i),tmpf(i); end;end; 下面編寫四個(gè)函數(shù)：function u = DecisF1( k,x ) %決策函數(shù)if k=4 u=x;else u=0:x;endfunction y = TransF1( k,x,u ) %狀態(tài)轉(zhuǎn)移方程y=x-u;function v = SubObjF1( k,x,u ) %階段k的指標(biāo)函數(shù)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)閙infunction y = ObjF1( v,f ) %基本方程中的函數(shù)y=v+f;y=-y; %max變?yōu)閙in測試代碼：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&