建立动态规划模型如下
代码使用matlab进行实现
代码部分先是编程了一个关于自由始段和终端的动态规划，求指标函数最小值的逆序算法递归计算程序：
%dynprog.m
function [p_opt,fval]=dynprog(x,DecisFun,SubObjFun,TransFun,ObjFun)
%DecisFun(k,x):由阶段k的状态变量x算出相应的允许决策变量
%SubObjFun(k,x,u):阶段指标函数
k=length(x(1,:)); %k的大小为x的整个第一行的长度即为x的列数
x_isnan=~isnan(x); %返回逻辑值
t_vub=inf; %创建所有值均为Inf的数组，inf表示正无穷大
t_vubm=inf*ones(size(x)); %创建与x相同尺寸的无穷大数组
f_opt=nan*ones(size(x)); %nan:非数值
d_opt=f_opt;
%计算终端（最后一阶段）相关值
tmp1=find(x_isnan(:,k)); %查找x_isnan整个第k列的非零元素
temp2=length(tmp1);
for i=1:temp2
u=feval(DecisFun,k,x(i,k));
tmp3=length(u);
for j=1:tmp3
tmp=feval(SubObjFun,k,x(tmp1(i),k),u(j));
if tmp<=t_vub
f_opt(i,k)=tmp;
d_opt(i,k)=u(j);
t_vub=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(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;
tmp60=find(tmp50==0);
if ~isempty(tmp60)
if nargin<5,tmp00=tmp00+f_opt(tmp60(1),ii+1);
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(tmp1,1);
fval=fval(find(~isnan(fval)),1);
%记录最优决策、最优轨线和相应指标函数值
p_opt=[];
tmpx=[];
tmpd=[];
tmpf=[];
tmp0=find(x_isnan(:,1));
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-1,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
再按题目编写对应程序：
主程序
%Q1.m
clear;
x=nan*ones(110,104); %x的最大范围[0,109]，因此x=0,1,...,109,所以x初始化取110行，列数表示共104周，nan表示无意义元素
for i=1:104
x(:,i)=(0:109)'; %按周定义x的可能取值
end
[p,f]=dynprog(x,@f1,@f2,@f3)
子程序
%f1.m
function u=f1(k,x) %在阶段k由状态变量x的值求出其相应的决策变量所有的取值
q=reshape(load('附件2.txt')',[1,104]);
if q(k)-x<0
u=0:109; %决策变量不能取为负值
else
u=q(k)-x:109; %产量满足需求且不超过109
end
u=u(:);
end
子程序
%f2.m
function v=f2(k,x,u) %阶段k的指标函数
q=reshape(load('附件2.txt')',[1,104]);
v=640*u+30*x+5*ceil(0.4*u)-30*q(k)+380;
end
子程序
%f3.m
function y=f3(k,x,u) %状态转移方程
y=x+u;
end
最后运行结果很奇怪，如图
p的值太多放不全
