动态规划在火力分配中的应用。

1. 问题描述
设有m 个目标，目标价值（重要性和危害性）各不相同，用数值A K （K=1，2，..m ）表示，计划用n 枚导弹突袭，导弹击毁目标的概率P K =1−e −a k u k ，其中a k 是常数，取决于导弹的特性与目标的性质；u k 为向目标发射的导弹数，问题：做出方案使预期的突击效果最大。

2. 问题建模
上述问题可以表述为
max V =∑A k (1−e
−a k u k )m k=1 约束条件为
∑u k m k=1=n （u k 为非负整数）
3. 算法描述
下面通过一个实例说明：设目标数目为4（m=4），导弹为5（n=5），A k 和a K 取值情况如下表所示：
表1：A k ，u k 取值情况
将火力分配可分为4个阶段，每个阶段指标函数为： V 1(u 1)=8(1−e −0.2u 1) V 2(u 2)=7(1−e −0.3u 2) V 3(u 3)=6(1−e −0.5u 3) V 4(u 4)=3(1−e −0.9u 4) u k 可能取值为0，1，2，3，4，5，将函数值带人如下表：
表2 函数值
动态规划问题基本方程为：
f k (x
k
)=max {V
k
(u
k
)+f
k+1
(x
k
−u k)}c
f 5(x
5
)=0
逐次向前推一级
K=4 f
4(x
4
)=V
4
(u
4
)=3(1−e−0.9u4)
K=3 f
3(x
3
)=max {V
3
(u
3
)+f
4
(x
3
−u3)}=max {6(1−e−0.5u3)+
f
4(x
3
−u3)}
K=2 f
2(x
2
)=max {V
2
(u
2
)+f
3
(x
2
−u2)}= max {7(1−e−0.3u2)+
f
3(x
2
−u2)}
K=1 f
1(x
1
)= max {V
1
(u
1
)+f
2
(x
1
−u1)}= max {8(1−e−0.2u1)+
f
2(x
1
−u1)}（0<u k<5可取等号）
只需要求解f
1
(5)的最大值然后反推回去就可以获得最优的分配方案
4.Matlab仿真求解
因为x k与u k取值为整数，可以采用动态规划的方法，获得f
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); % 最后一步状态个数
for i=1:tmp2
u=feval(DecisFun,k,x(tmp1(i),k));
tmp3=length(u);%决策变量
for j=1:tmp3 % 求出当前状态下所有决策的最小性能指标
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 %求出当前状态下所有可能的决策
u=feval(DecisFun,ii,x(tmp10(i),ii));
tmp30=length(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) ;
if~isempty(tmp60)
if nargin<6 %矩阵不同需要修改nargin的值，很重要
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 %当前状态的性能指标
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;
end
fval=f_opt(:,1);
tmp0 = find(~isnan(fval));
fval=fval(tmp0,1);
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 ~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;
下面编写四个函数：
function u = DecisF1( k,x ) %决策函数
if k==4
u=x;
else
u=0:x;
end
function y = TransF1( k,x,u ) %状态转移方程
y=x-u;
function v = 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变为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','SubObjF1','TransF1','ObjF1')
%p为分配方案，f为结果
5.Matlab仿真结果分析
运行结果显示：
P为方案：
即对目标1发射1枚导弹，对目标2发射1枚，对目标3发射2枚，对目标4发生1枚导弹，能获得最大效能。

f为最大的效能：
即在这种分配下的最大效能为8.8374。