数学建模作业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

%% 逆推计算各阶段的递归调用程序

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

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

function u=DecisFun(k,x)