动态规划matlab仿真实例

动态规划在火力分配中的应用。1.问题描述设有m个目标，目标价值（重要性和危害性）各不相同，用数值AK（K=1，2，..m）表示，计划用n枚导弹突袭，导弹击毁目标的概率PK=1−𝑒−𝑎𝑘𝑢𝑘，其中𝑎𝑘是常数，取决于导弹的特性与目标的性质；𝑢𝑘为向目标发射的导弹数，问题：做出方案使预期的突击效果最大。2.问题建模上述问题可以表述为max𝑉=∑𝐴𝑘(1−𝑒−𝑎𝑘𝑢𝑘)𝑚𝑘=1约束条件为∑𝑢𝑘𝑚𝑘=1=𝑛（𝑢𝑘为非负整数）3.算法描述下面通过一个实例说明：设目标数目为4（m=4），导弹为5（n=5），𝐴𝑘和aK取值情况如下表所示：表1：Ak，𝑢𝑘取值情况目标K1234𝐴𝑘8763𝑢𝑘0.20.30.50.9将火力分配可分为4个阶段，每个阶段指标函数为：𝑉1(𝑢1)=8(1−𝑒−0.2𝑢1)𝑉2(𝑢2)=7(1−𝑒−0.3𝑢2)𝑉3(𝑢3)=6(1−𝑒−0.5𝑢3)𝑉4(𝑢4)=3(1−𝑒−0.9𝑢4)𝑢𝑘可能取值为0，1，2，3，4，5，将函数值带人如下表：表2函数值u𝑉1(𝑢1)𝑉2(𝑢2)𝑉3(𝑢3)𝑉4(𝑢4)0000011.451.812.361.7922.643.163.792.5133.614.154.662.8144.414.895.192.9355.065.445.512.97动态规划问题基本方程为：𝑓𝑘(𝑥𝑘)=max{𝑉𝑘(𝑢𝑘)+𝑓𝑘+1(𝑥𝑘−𝑢𝑘)}c𝑓5(𝑥5)=0逐次向前推一级K=4𝑓4(𝑥4)=𝑉4(𝑢4)=3(1−𝑒−0.9𝑢4)K=3𝑓3(𝑥3)=max{𝑉3(𝑢3)+𝑓4(𝑥3−𝑢3)}=max{6(1−𝑒−0.5𝑢3)+𝑓4(𝑥3−𝑢3)}K=2𝑓2(𝑥2)=max{𝑉2(𝑢2)+𝑓3(𝑥2−𝑢2)}=max{7(1−𝑒−0.3𝑢2)+𝑓3(𝑥2−𝑢2)}K=1𝑓1(𝑥1)=max{𝑉1(𝑢1)+𝑓2(𝑥1−𝑢1)}=max{8(1−𝑒−0.2𝑢1)+𝑓2(𝑥1−𝑢1)}（0𝑢𝑘5可取等号）只需要求解𝑓1(5)的最大值然后反推回去就可以获得最优的分配方案4.Matlab仿真求解因为𝑥𝑘与𝑢𝑘取值为整数，可以采用动态规划的方法，获得𝑓1(5)的最大值，对应的最优方案function[p_opt,fval]=dynprog(x,DecisFun,SubObjFun,TransFun,ObjFun)%求解动态规划问题最小值函数k=length(x(1,:))%判断决策级数x_isnan=~isnan(x);%非空状态矩阵t_vubm=inf*ones(size(x));%性能指标中间矩阵f_opt=nan*ones(size(x));%总性能指标矩阵d_opt=f_opt;%每步决策矩阵tmp1=find(x_isnan(:,k));%最后一步状态向量tmp2=length(tmp1);%最后一步状态个数fori=1:tmp2u=feval(DecisFun,k,x(tmp1(i),k));tmp3=length(u);%决策变量forj=1:tmp3%求出当前状态下所有决策的最小性能指标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%求出当前状态下所有可能的决策u=feval(DecisFun,ii,x(tmp10(i),ii));tmp30=length(u);forj=1:tmp30%求出当前状态下所有决策的最小性能指标tmp00=feval(SubObjFun,ii,x(tmp10(i),ii),u(j));%单步性能指标tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j));%下一状态tmp50=x(:,ii+1)-tmp40;%找出下一状态在x矩阵的位置tmp60=find(tmp50==0);if~isempty(tmp60)ifnargin6%矩阵不同需要修改nargin的值，很重要tmp00=tmp00+f_opt(tmp60(1),ii+1);%setthedefaultobjectvalueelsetmp00=feval(ObjFun,tmp00,f_opt(tmp60(1),ii+1));end%当前状态的性能指标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(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);if~isempty(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(i),tmpd(i),tmpf(i)];end;end;下面编写四个函数：functionu=DecisF1(k,x)%决策函数ifk==4u=x;elseu=0:x;endfunctiony=TransF1(k,x,u)%状态转移方程y=x-u;functionv=SubObjF1(k,x,u)%阶段k的指标函数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变为minfunctiony=ObjF1(v,f)%基本方程中的函数y=v+f;y=-y;%max变为min测试代码：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为结果5.Matlab仿真结果分析运行结果显示：P为方案：即对目标1发射1枚导弹，对目标2发射1枚，对目标3发射2枚，对目标4发生1枚导弹，能获得最大效能。f为最大的效能：即在这种分配下的最大效能为8.8374