现在时间是:
当前位置:首 页 >> 数学建模>> 教学区>> 文章列表

数学建模作业10

作者:芮逸文   发布时间:2020-05-13 22:26:50   浏览次数:106

 

1、最短线路问题

clear

clc

% user input data

% 'a':每个阶段的站点数目

% 'b':所有站点的距离

% num:所有站点的个数总和

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;

 

% c:最小路径

% c1:第一个决策点所产生的四个距离

% c2:其他决策点所产生的四个距离,并将较小值放入c1

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、生产——存贮问题

% eg13f1_2.m

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

 

% eg13f2_2.m

 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

 

 

% eg13f3_2.m

 function s_next=TransFun(k,x,u)

d = [ 2 3 2 4 ];

s_next = x + u - d(k);

 

% DynProg.m

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

 

%main.m

clear all

clc

 

x1 = 0:4;

s=nan*ones(5,1);

s(1) = 0;

x = [s x1' x1' x1'];

[p_opt,fval]=dynprog(x,@eg13f1_2,@eg13f2_2,@eg13f3_2);

fval

fval

 

运行结果:

fval =

 

   20.5000

 

 

 

 

3、生产——库存管理系统

% eg13f1_2.m

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

 

% eg13f2_2.m

 function V=ObjFun(k,x,u)

V = x+ 0.005*u^2;

 

% eg13f3_2.m

 function s_next=TransFun(k,x,u)

d = [ 600 700 500 1200];

s_next = x + u - d(k);

 

%dynprog.m

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

 

%main.m

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,@eg13f1_2,@eg13f2_2,@eg13f3_2);

fval

 

运行结果:

fval =

 

       11800

 

 

4、企业生产管理问题

% eg13f1_2.m

function u=DecisFun(k,x)

u = 0 : x;

 

% eg13f2_2.m

 function V=ObjFun(k,x,u)

if u>900

        V = 0;

else

        V =3*u+5*x;

end

 

% eg13f3_2.m

function s_next=TransFun(k,x,u)

a = 0.7;

b = 0.9;

s_next =(a - b)*u + b*x;

 

%dynprog.m

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

 

%main.m

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,@eg13f1_2,@eg13f2_2,@eg13f3_2);

fval

 

运行结果:

fval =

 

       12052

 

5、设备更新问题

% eg13f1_2.m

function u=DecisFun(k,x)

u = 1:x;

 

% eg13f2_2.m

 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 = 5-0.5 - C(x);

end

 

% eg13f3_2.m

function s_next=TransFun(k,x,u)

if u ==6

    s_next = x + 1;

else

    s_next = 1;

end

 

%dynprog.m

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

 

%main.m

clear all

clc

x1 = (1:5);

s=nan*ones(5,1);

s(1) = 5;

x = [s x1' x1' x1' x1'];

[p_opt,fval]=dynprog(x,@eg13f1_2,@eg13f2_2,@eg13f3_2);

fval

 

运行结果:

fval =

 

   18.5000







上一篇:没有了    下一篇:没有了

Copyright ©2020    计算数学达人 All Right Reserved.

技术支持:自助建站 | 领地网站建设 |短信接口 |燕窝 版权所有 © 2005-2020 lingw.net.粤ICP备16125321号 -5