数学建模第十次作业
作者：杨创 发布时间：2020-05-14 09:11:07 浏览次数：97例1（最短路线问题）
程序：
a=[1 2 4 3 3 2 1];
b=[5 3 -1 -1;
1 3 6 -1;-1 8 7 6;
6 8 -1 -1;3 5 -1 -1;-1 3 3 -1;-1 8 4 -1;
2 2 -1 -1;-1 1 2 -1;-1 3 3 -1;
3 5 -1 -1;5 2 -1 -1;6 6 -1 -1;
4 -1 -1 -1;3 -1 -1 -1];
num=16;
max_col=4;
a_length=size(a);
route=zeros(a_length(2),max_col);
row=2;
c=b(1,:);
c1=zeros(1,max_col);c2=zeros(1,max_col);
for i=2:a_length(2)-1
n=a(i);
for k=1:max_col
if b(row,k)>0
c1(k)=b(row,k)+c(1);
route(i,k)=1;
end
end
if n>1
for j=1:n-1
for k=1:max_col
if b(row+j,k)>0
c2(k)=b(row+j,k)+c(j+1);
if c2(k)<c1(k)||c1(k)==0
c1(k)=c2(k);
route(i,k)=j+1;
end
end
end
end
end
row=row+n;
c=c1;
c1=zeros(1,max_col);c2=zeros(1,max_col);
end
correct_col=route(end-1,1);
l_route=route(:,correct_col);
l_route(1)=1;
l_route(end-1)=correct_col;
l_route(end)=1;
disp('最小路径为');disp(c(1));disp('路径规划为');disp(l_route');
结果：
最小路径为
18
路径规划为
1 1 2 1 2 2 1
例2（生产—存贮问题）
程序：
function [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)
k=length(x(1,:)); % 判读决策级数
f_opt=nan*ones(size(x)); % 总性能指标
d_opt=f_opt; % 每步决策
t_vubm=inf*ones(size(x)); % 性能指标，中间矩阵
x_isnan=~isnan(x); % 非空状态矩阵
% t_vub=inf;
%% 计算终端相关值
tmp1=find(x_isnan(:,k)); % 最后一步状态状态向量
tmp2=length(tmp1); % 最后一步状态状态个数
for i=1:tmp2
u=DecisFun(k,x(i,k));
tmp3=length(u);
for j=1:tmp3
tmp=ObjFun(k,x(tmp1(i),k),u(j));
if tmp<=t_vubm(i,k)
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=DecisFun(ii,x(i,ii));
tmp30=length(u);
for j=1:tmp30 % 求出当前状态下所有决策的最小性能值
tmp00=ObjFun(ii,x(tmp10(i),ii),u(j));% 单步性能指标
tmp40=TransFun(ii,x(tmp10(i),ii),u(j)); % 下一状态
tmp50=x(:,ii+1)-tmp40; % 找出下一状态在矩阵中的位置
tmp60=find(tmp50==0);
if ~isempty(tmp60)
tmp00=tmp00+f_opt(tmp60(1),ii+1);
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);
%% 记录最优决策、最优轨线和相应指标函数值
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)=ObjFun(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)=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)=ObjFun(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 V=ObjFun(k,x,u)
d = [ 2 3 2 4 ];
if u == 0
V = 0.5 * (x + u -d(k));
else
if u>6
V = 10^6;
else
V = 3 + u +0.5*(x + u -d(k));
end
end
function s_next=TransFun(k,x,u)
d = [ 2 3 2 4 ];
s_next = x + u - d(k);
function u=DecisFun(k,x)
m=6;
d = [ 2 3 2 4 ];
if k == 4
u = d(k) - x;
else
u = max (0,d(k)-x):m;
end
x1 = 0:4;
s=nan*ones(5,1);
s(1) = 0;
x = [s x1' x1' x1'];
[p_opt,fval]=dynprog(x,@DecisFun,@ObjFun,@TransFun)
结果：
p_opt =
1.0000 0 5.0000 9.5000
2.0000 3.0000 0 0
3.0000 0 6.0000 11.0000
4.0000 4.0000 0 0
fval =
20.5000
NaN
NaN
NaN
NaN
例3（生产—库存管理系统）
程序：
function u=DecisFun(k,x)
d =[ 600 700 500 1200];
if k == 4
u = d(k) - x;
else
u = max (0,d(k)-x):1200;
end
function V=ObjFun(k,x,u)
V = x+ 0.005*u^2;
function s_next=TransFun(k,x,u)
d = [ 600 700 500 1200];
s_next = x + u - d(k);
x1 =100*(0:4);
s = nan*ones(5,1);
s(1)=0;
x = [s x1' x1' x1'];
[p_opt,fval]=dynprog(x,@DecisFun,@ObjFun,@TransFun)
结果：
p_opt =
1 0 600 1800
2 0 700 2450
3 0 800 3200
4 300 900 4350
fval =
11800
NaN
NaN
NaN
NaN
例4（企业生产管理问题）
程序：
function V=ObjFun(k,x,u)
V =8*u + 5*(x-u);
End
function s_next=TransFun(k,x,u)
a = 0.7;
b = 0.9;
s_next =(a - b)*u + b*x;
function u=DecisFun(k,x)
u = 0 : x;
x1 = (500:999);
s=nan*ones(500,1);
s(1) = 1000;
x = [s x1' x1' x1' x1'];
[p_opt,fval]=dynprog(x,@DecisFun,@ObjFun,@TransFun)
结果：
p_opt =
1 1000 1000 8000
2 700 10 3530
3 628 6 3158
4 564 3 2829
5 507 0 2535
fval =
20052
例5（设备更新问题）
function V=ObjFun(k,x,u)
R = [ 5 4.5 4 3.75 3 2.5]; % 产生的经济效益
U = [ 0.5 1 1.5 2 2.5 3]; % 维修费用
C = [ 0.5 1 1.5 2 2.5 3]; % 换购一台新的设备所需要的费用
if u == 1
V = R(x) - U(x);
else
V = R(1) - U(1) - C(x);
End
function s_next=TransFun(k,x,u)
if u == 1
s_next = x + 1;
else
s_next = 1;
end
function u=DecisFun(k,x)
u = [1 0];
x1 = (1:6);
s=nan*ones(6,1);
s(1) = 1;
x = [s x1' x1' x1' x1'];
[p_opt,fval]=dynprog(x,@DecisFun,@ObjFun,@TransFun)
结果：
p_opt =
1.0000 1.0000 1.0000 4.5000
2.0000 2.0000 1.0000 3.5000
3.0000 3.0000 1.0000 2.5000
4.0000 4.0000 1.0000 1.7500
5.0000 5.0000 1.0000 0.5000
fval =
12.7500
NaN
NaN
NaN
NaN
NaN
