数学建模作业10
作者:18074017闫培韬 发布时间:2020-05-14 12:13:39 浏览次数:90
1.function u=DecisFun(k,x)
if x == 1 ,u = [2;3];
elseif x == 2,u = [4;5;6];
elseif x == 3,u = [5;6;7];
elseif (x==4)||(x==5), u = [8;9];
elseif (x==6)||(x==7),u = [9;10];
elseif x==8,u = [11,12];
elseif (x==9)||(x==10),u = [12;13];
elseif (x==11)||(x==12)||(x==13),u = [14;15];
elseif (x==14)||(x==15),u = 16;
elseif x==16,u =16;
end
clear all
clc
x=nan*ones(4,7);
x(1,1)=1;
x(1:2,2) = [2;3];
x(1:4,3)=(4:7)';
x(1:3,4)=(8:10)';
x(1:3,5)=(11:13)';
x(1:2,6)=[14;15];
x(1,7)=16;
[p_opt,fval]=dynprog(x,@DecisFun,@ObjFun,@TransFun);
function V=ObjFun(k,x,u)
tt=[5;3;1;3;6;8;7;6;6;8;3;5;3;3;8;4;2;2;1;2;3;3;3;5;5;2;6;6;4;3];
tmp = [x == 1&& u == 2,x == 1&& u == 3,x == 2&& u == 4,x == 2&& u == 5,...
x == 2&& u == 6,x == 3&& u == 5,x == 3&& u == 6,x == 3&& u == 7,...
x == 4&& u == 8,x == 4&& u == 9,x == 5&& u == 8,x == 5&& u == 9,...
x == 6&& u == 9,x == 6&& u == 10,x == 7&& u == 9,x == 7&& u == 10,...
x == 8&& u == 11,x == 8&& u == 12,x == 9&& u == 11,x == 9&& u == 12,...
x == 10&& u == 12,x == 10&& u == 13,x == 11&& u == 14,x == 11&& u == 15,...
x == 12&& u == 14,x == 12&& u == 15,x == 13&& u == 14,x == 13&& u == 15,...
x == 14&& u == 16,x == 15&& u == 16];
V = tmp*tt;
function s_next=TransFun(k,x,u)
s_next=u;
2.function u=DecisFun(k,x)
d =100* [ 6 7 5 12];
if k == 4
u = d(k) - x;
else
u = max (0,d(k)-x):1200;
end
clear all
clc
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);
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);
function [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)
% [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)
% 自由始端和终端的动态规划,求指标函数最小值的逆序算法递归
% 计算程序。x是状态变量,一列代表一个阶段状态;M-函数
% DecisFun(k,x)由阶段k的状态变量x求出相应的允许决策变量;
% M-函数ObjFun(k,x,u)是阶段指标函数,M-函数TransFun(k,x,u)
% 是状态转移函数,其中x是阶段k的某状态变量,u是相应的决策变量;
% 输出p_opt由4列构成,p_opt=[序号组;最优策略组;最优轨线组;
% 指标函数值组];fval是一个列向量,各元素分别表示p_opt各
% 最优策略组对应始端状态x的最优函数值;
%先写3个函数
% eg13f1_2.m
% function u=DecisF_1(k,x)
% 在阶段k由状态变量x的值求出其相应的决策变量所有的取值
% c=[70,72,80,76];q=10*[6,7,12,6];
% if q(k)-x<0,u=0:100; %决策变量不能取为负值
% else,u=q(k)-x:100;end; %产量满足需求且不超过100
% u=u(:);
% eg13f2_2.m
% function v=ObjF_1(k,x,u)
% 阶段k的指标函数
% c=[70,72,80,76];v=c(k)*u+2*x;
% eg13f3_2.m
% function y=TransF_1(k,x,u)
% 状态转移方程
% q=10*[6,7,12,6];y=x+u-q(k);
%调用DynProg.m计算如下:
% clear;x=nan*ones(14,4);% x是10的倍数,最大范围0≤x≤130,
% %因此x=0,1,...13,所以x初始化取14行,nan表示无意义元素
% x(1:7,1)=10*(0:6)'; % 按月定义x的可能取值
% x(1:11,2)=10*(0:10)';x(1:12,3)=10*(2:13)';
% x(1:7,4)=10*(0:6)';
% [p,f]=dynprog(x,'eg13f1_2','eg13f2_2','eg13f3_2')
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 u=DecisFun(k,x)
if x == 1 ,u = [2;3];
elseif x == 2,u = [4;5;6];
elseif x == 3,u = [5;6;7];
elseif (x==4)||(x==5), u = [8;9];
elseif (x==6)||(x==7),u = [9;10];
elseif x==8,u = [11,12];
elseif (x==9)||(x==10),u = [12;13];
elseif (x==11)||(x==12)||(x==13),u = [14;15];
elseif (x==14)||(x==15),u = 16;
elseif x==16,u =16;
end
function V=ObjFun(k,x,u)
tt=[5;3;1;3;6;8;7;6;6;8;3;5;3;3;8;4;2;2;1;2;3;3;3;5;5;2;6;6;4;3];
tmp = [x == 1&& u == 2,x == 1&& u == 3,x == 2&& u == 4,x == 2&& u == 5,...
x == 2&& u == 6,x == 3&& u == 5,x == 3&& u == 6,x == 3&& u == 7,...
x == 4&& u == 8,x == 4&& u == 9,x == 5&& u == 8,x == 5&& u == 9,...
x == 6&& u == 9,x == 6&& u == 10,x == 7&& u == 9,x == 7&& u == 10,...
x == 8&& u == 11,x == 8&& u == 12,x == 9&& u == 11,x == 9&& u == 12,...
x == 10&& u == 12,x == 10&& u == 13,x == 11&& u == 14,x == 11&& u == 15,...
x == 12&& u == 14,x == 12&& u == 15,x == 13&& u == 14,x == 13&& u == 15,...
x == 14&& u == 16,x == 15&& u == 16];
V = tmp*tt;
function s_next=TransFun(k,x,u)
s_next=u;
clear all
clc
x=nan*ones(4,7);
x(1,1)=1;
x(1:2,2) = [2;3];
x(1:4,3)=(4:7)';
x(1:3,4)=(8:10)';
x(1:3,5)=(11:13)';
x(1:2,6)=[14;15];
x(1,7)=16;
[p_opt,fval]=dynprog(x,@DecisFun,@ObjFun,@TransFun);
3. function u=DecisFun(k,x)
d =100* [ 6 7 5 12];
if k == 4
u = d(k) - x;
else
u = max (0,d(k)-x):1200;
end
clear all
clc
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);
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);
4. function u=DecisFun(k,x)
u = 0 : x;
clear all
clc
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);
function V=ObjFun(k,x,u)
if u>900
V = 0;
else
V =8*u + 5*(x-u);
V=-V;
End
function s_next=TransFun(k,x,u)
a = 0.7;
b = 0.9;
s_next =(a - b)*u + b*x;
5. function u=DecisFun(k,x)
u = [1 0];
clear all
clc
x1 = (1:5);
s=nan*ones(5,1);
s(1) = 1;
x = [s x1' x1' x1' x1'];
[p_opt,fval]=dynprog(x,@DecisFun,@ObjFun,@TransFun);
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 = 5-[ 4.5 4 3.5 3 2.5 2]; % 换购一台新的设备所需要的费用
if u == 1
V = R(x) - U(x);
else
V = R(x) - U(x) - C(x);
End
function s_next=TransFun(k,x,u)
if u == 1
s_next = x + 1;
else
s_next = 1;
end