动态规划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 ，𝑢𝑘取值情况

目标K 𝐴𝑘 𝑢𝑘 1 8 2 7 3 6 4 3 将火力分配可分为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) 0 1 2 3 4 5 0 0 0 0 动态规划问题基本方程为：

𝑓𝑘(𝑥𝑘)=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); % 最后一步状态个数 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=[,,,]; 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为最大的效能：

即在这种分配下的最大效能为