动态规划matlab仿真实例

1、动态规划在火力分配中的应用。1.问题描述设有m个目标，目标价值(重要性和危害性)各不相同，用数值A(K=1, 2, .m) 表示，计划用n枚导弹突袭，导弹击毁目标的概率Pk=1- ?- ?，其中？?是 常数，取决于导弹的特性与目标的性质；？?为向目标发射的导弹数，问题：做 出方案使预期的突击效果最大。2问题建模上述问题可以表述为max? = E?=1?(1- ?- ?)约束条件为?1 ? = ? （?为非负整数）3. 算法描述下面通过一个实例说明：设目标数目为 4 (m=4,导弹为5 (n=5)，?和aK取 值情况如下表所示：表1: A，?取值情况目标K1234?8763?将火力分配可分为4个

2、阶段，每个阶段指标函数为：?1(?1)=8(1 -?- 0.2?1)?2(?2)= 7(1 -?- 3?2)?3(?3)=6(1 -?- 0.5?3)?4(?4)= 3(1 -?- 9?4)?可能取值为0, 1, 2，3, 4, 5,将函数值带人如下表：表2函数值u?1（?1）?2（ ?2）?3（?3）?4（ ?4）0000012345动态规划问题基本方程为：?(?) = ma?(?) + ?+i(?- ?) c?5(?5)=0逐次向前推一级K=4?4( ?4)= ?4( ?4)= 3(1 - ?-0.9?4)K=3?3(?3)= ma?3(?3)+ ?4(?3 - ?3) = ma?6(1

3、-?一05?3)+?4(?3 - ?3)K=2?2(?2)= ma?2(?2)+ ?3(?2 - ?2) = ma?7(1 -?-0.3?Q)+?3(?2 - ?2)K=1?i(?i) = maX?i(?i) + ?2(?1 - ?i) = maW8(1 -?-02?1) +?2(?1 - ?i)(0 ? 5可取等号)只需要求解？?i(5)的最大值然后反推回去就可以获得最优的分配方案4. Matlab仿真求解因为？?？与？?取值为整数，可以采用动态规划的方法，获得？?i(5)的最大值,对应的最优方案求解动态规fun cti on p_opt,fval=d yn prog(x,DecisFu n

4、, SubObjFu n,Tran sFu n,ObjFun) %划问题最小值函数k=le ngth(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); % 最后一步状态个数for i=1:tmp2u=feval(DecisFun,k,x(tmp1(i),k);tmp3=length(u);%

5、 决策变量for j=1:tmp3 % 求出当前状态下所有决策的最小性能指标tmp=feval(SubObjFun,k,x(tmp1(i),k),u(j);if tmp = t_vubm(i,k) %t_vubf_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vubm(i,k)=tmp;end;end;endfor ii=k-1:-1:1tmp10=find(x_isnan(:,ii);tmp20=length(tmp10);for i=1:tmp20 % 求出当前状态下所有可能的决策u=feval(DecisFun,ii,x(tmp10(i),ii);tmp30=length

6、(u) ;for j=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) ;ifisempty(tmp60)if nargin6 % 矩阵不同需要修改 nargin 的值，很重要 tmp00=tmp00+f_opt(tmp60(1),ii+1); % set the d

7、efault object valueelsetmp00=feval(ObjFun,tmp00,f_opt(tmp60(1),ii+1);end % 当前状态的性能指标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); p_opt=;tmpx=;tmpd=;tmpf=; tmp01=length(tmp0);for i=1:tmp01tmpd(

8、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 isempty(tmp2)tmpd(i)=d_opt(tmp2(1),ii);endtmpf(i)=feval(SubObjFun,i

9、i,tmpx(i),tmpd(i);p_opt(k*(i-1)+ii,1,2,3,4)=ii,tmpx(i),tmpd(i),tmpf(i);end;end;决策函数状态转移方程阶段k的指标函数基本方程中的函数下面编写四个函数： function u = DecisF1( k,x ) % if k=4u=x;elseu=0:x;endfunction y = TransF1( k,x,u ) % y=x-u;function v = SubObjF1( k,x,u ) % a=,;A=8,7,6,3; v=A(k)*(1-exp(-a(k)*u); v=-v; %max 变为 min function y = 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,