opt_obj=cur_obj;
opt_dec=x;
end
end%forx4的循环结束
end%forx3的循环结束
end%forx2的循环结束
end%forx1的循环结束
disp（sprintf（'最优目标＝%20.6f',opt_obj））
fori=1:
N,
disp（sprintf（'x（%d）=%20.6f',i,x（i）））
end
exp_time=etime（clock,time_begin）
C=ones（3,1）;%相对于讲目标函数中单价设为1，则计算结果应为总长
度S1+S2+S3+S4+S5
len=objfun_97（opt_dec）
程序运行结果:
最优目标＝22603.376739
x
（1）=30.000000；x
（2）=30.000000；x（3）=30.000000；
x（4）=30.000000。
exp_time=105.7900
len=39.3180
案例：
背包问题
有一组物品S，共有9件，其中第i件重
，价值
，从S中取出一些物品出来装背包，使总价值最大，而不超过总重量的给定上限30kg。
i
1
2
3
4
5
6
7
8
9
（kg）
2
1
1
2.5
10
6
5
4
3
（元）
10
45
30
100
150
90
200
180
300
问题分析
这是一个典型的最优化问题，优化目标是总价值最大，决策是决定装哪些物品，而装载物品又受到背包所能承受重量30kg的限制。
因此可以建立该问题的最优化数学模型，而且是0－1整数规划模型。
变量与符号说明
：
用来表示是否装载第i件物品，如果
表示不装载该物品，如果
表示装载该物品（
），令
。
：
第i件物品重量（
），单位：
kg。
：
第i件物品价值（
），单位：
元。
：
表示按照决策
，得到的装载物品的总价值，单位：
元。
模型建立
通过前面的分析及变量的定义，可得
。
原问题抽象为如下的0－1整数线性规划模型：
模型的求解
（1）Lingo程序：
MODEL:
!
求解背包;
!
定义所需集合和数组变量;
SETS:
SERIES/1..9/:
weight,value,xx;
ENDSETS
!
输入题上所给数据;
DATA:
weight=2,1,1,2.5,10,6,5,4,3;
value=10,45,30,100,150,90,200,180,300;
ENDDATA
MAX=@SUM（SERIES（i）:
value（i）*xx（i））;
@Sum（SERIES（i）:
weight（i）*xx（i））<=30;
@for（SERIES（i）:
@bin（xx（i）））;
END
结果为：
最大价值为：
1015.000，总质量=28.50kg取其中的第1、2、3、4、5、7、8、9件物品放入背包中。
（2）贪婪算法用Matlab语言求解：
w=[2112.5106543];
v=[10453010015090200180300];
g=v./w;%单位质量的价值（指标）
maxw=30;
%对指标排序，以便根据指标选择物品
[y,idx]=sort（g）;
curw=0;%存储当前装的总重
curv=0;%存储当前装的价值
hasidx=[];%存储当前装的物品序号
fori=length（g）:
-1:
1,%从指标g最大的开始选择
id=idx（i）;%important
ifcurw+w（id）<=maxw,
curw=curw+w（id）;
curv=curv+v（id）;
hasidx=[hasidx,id];
else
break;%结束，不能再装
end
end
disp（sprintf（'总共装载物品数=%5d件',length（hasidx）））
disp（sprintf（'质量=%10.2fkg',curw））
disp（sprintf（'价值=%10.2f元',sum（v（hasidx））））
结果：
最大总价值为945.00元，总质量=22.50kg。
其中第9、8、2、7、4、3、6件物品放入背包中。
（3）穷举法用Matlab求解：
%背包问题求解得贪婪算法实现
w=[2112.5106543];%各件物品重量
v=[10453010015090200180300];%价值
maxw=30;%背包承受重量上限
optvalue=-1;%最优目标（最大价值）初始化
forx1=0:
1,
forx2=0:
1,
forx3=0:
1,
forx4=0:
1,
forx5=0:
1,
forx6=0:
1,
forx7=0:
1,
forx8=0:
1,
forx9=0:
1,
x=[x1x2x3x4x5x6x7x8x9];%决策向量
%不超过重量限制，也要比以前装包方案好
if（dot（x,w）<=maxw）&（dot（x,v）>optvalue）,
optx=x;%存储
optvalue=dot（x,v）;%存储
end
end
end
end
end
end
end
end
end
end
hasidx=find（optx~=0）;
optx
disp（sprintf（'总共装载物品数=%5d件',length（hasidx）））
disp（sprintf（'质量=%10.2fkg',dot（optx,w）））
disp（sprintf（'价值=%10.2f元',dot（optx,v）））
案例：
运载问题
有一个最大货运量为30t的卡车，用来装载三种货物。
每种货物的单位重量及相应的单位价值如表所示。
问应如何装载才能使总的价值最大？
货物编号
1
2
3
单位重量
3
4
5
单位价值
4
5
6
问题分析
这是一个典型的最优化问题，优化目标是卡车装载的总价值最大，装载当然越多越好。
但又受到卡车最大载货量的限制。
因此可以建立该问题的最优化数学模型，而且是整数规划模型。
变量与符号说明
（
）：
卡车装载第i种货物的件数，为整数。
（
）：
第i种货物的单位重量，单位：
吨。
（
）：
第i种货物的单位价值，单位：
千元。
（
）：
卡车装载前i种货物的货运量。
单位：
吨
：
表示按照决策
，得到的装载物品的总价值，单位：
元。
模型建立
（1）
（1） 线性规划模型
通过前面的分析及变量的定义，可得
。
原问题抽象为如下整数线性规划模型：
此可以用前面的线性规划和非线性规划方法来求解
（2）动态规划模型（前向算法）
设
为当背包中允许装入货物的总重量不超过
，并且只允许装入前
种货物，采用最优策略时的最大价值。
则此问题的递推方程为
初始条件为：
状态转移方程为：
Matlab求解程序为：
function[value]=MyFirstDP
%动态规划模型求解示范
globalM;%矩阵的第i列针对第i种货物
M=[345;%第i种货物的单位重量
456];%第i种货物的单位价值
value=getmaxvalue（3,30）;
被调的子函数
functionvalue=getmaxvalue（k,y）
%这是一个子函数用于递归求最优值
globalM;
value=-inf;
ifk==1
value=M（2,1）*fix（y/M（1,1））;
which=fix（y/M（1,1））;
disp（sprintf（'第%d阶段f%d（%8.2f）=%8.2f=%d（单）*%d（决策）',…
k,k,y,value,M（2,k）,which））
return
end
fornum=0:
fix（y/M（1,k））
curvalue=M（2,k）*num+getmaxvalue（k-1,y-M（1,k）*num）;
ifcurvalue>value
value=curvalue;
which=num;
end
end
disp（sprintf（'第%d阶段f%d（%8.2f）=%8.2f=%d（单）*%d（决策）+f（%d,%8.2f）',…
k,k,y,value,M（2,k）,which,k-1,y-M（1,k）*which））
运行结果为：
ans=40，决策变量为：
。
（注：
大家可以照些程序来设计后向处法）
∙连续系统实验
实验目的
学会对实际问题进行数学抽象。
掌握常用的机理分析方法。
学会用相关软件对连续系统进行精确符号运算和数值计算。
实验原理
常用方法与思想：
平衡与增长：
许多问题在数量上的变化可以归结为促进因素和抑制因素，进而用此关系来建立变量间的关系。
类比：
利用问题相同的或想类似的机理，把问题归结或转化为我们熟知的模型。
逻辑演绎：
利用的规律，对实际问题进行适当的简化和变换，通过逻辑推理的演绎，以其满足的规律来建立变量之间的关系。
一般步骤：
1）1） 分析问题、识别问题，找出问题的主要因素和次要因素，确定变量和基本规律的类型。
先以主要因素为重点考察对象。
2）2） 用假设将次要因素简化，可以用常识或问题分析需要来确定对原问题简化的程度℃。
3）3） 确定变量间的相互关系，常以考察变量所满足的方程或动力系统的方式表述变量间的因果关系。
4）4） 用演绎的方式或数值方法式来求解或解释模型。
5）5） 选择合适的语言，编制程序求解。
6）6） 维修模型，把结果与实际比较，用前面得到的次要因素及简化的假设逐步放宽或加入，得到进一步的演绎和数值解，直到模型满足实际问题的需要和精度。
案例：
冷却问题
一个较热的物体置于室温为18℃的房间内，该物体最初的温度是60℃，3分钟以后降到50℃。
想知道它的温度降到30℃需要多少时间？
10分钟以后它的温度是多少？
问题分析
这是一个带有许多不定因素问题。
首先物体的形状不知道，室温条件是否变化不知道，热在物体内部的分析不知道，热的传播有幅射、传导、对流三种不同的方式，等等。
我们建立的模型有可能是偏微分方程。
为简化问题，可以认为物体每一点的温度都一样，只考虑传导过程，室温在冷却过程中保持不变，热交换只在物体与空气的接触面进行，而且在接触面两侧的温度差就是物体与空气的温度差。
基本假设
（1）假设房间足够大，放入温度较低或较高的物体时，室内温度基本不受影响。
（2）物体各点的温度总是保持一致。
（3）只考虑热传导过程。
（4）以放置物体时刻为记时初始时刻，时间以分钟为单位。
变量说明
名称
变量符号
单位
时间
t
分
室内温度
m
℃
物体的温度
T（t）
℃
建立模型
我们已知，在物理学中有牛顿冷却（加热）定律：
将温度为T的物体放入处于常温m的介质中时，T的变化速率正比于T与周围介质的温度m之差。
所以建立微分方程
其中参数k>0，室温m=18。
求解
（1）
（1） 用Mathematica语言求解
（*初始化变量的数据*）
m=.;k=.;T=.;
（*求解微分方程*）
sol=DSolve[{T'[t]==-k（T[t]-m）,T[0]==60},T[t],t];
TT=T[t]/.sol[[1]];
（*代入室内温度值*）
TT=TT/.m->18;
TK=TT/.t->3;
（*确定冷却系数k*）
kk=FindRoot[TK-50==0,{k,0.1}];
k=k/.kk;
（*求10分钟后，物体的温度*）
td=TT/.t->10
（*求当物体的温度为30时，经过的时间*）
FindRoot[TT==30,{t,13}]
结果为：
10分钟后，物体温度为35℃左右。
当物体的温度为34℃时，经过的时间为13.8分左右。
（2）
（2） 用Matlab语言求解
%初始化变量的数据
symsTtt0km;
%求解微分方程
yy=dsolve（'DT=-k*（T-m）','T（0）=60','t'）;
%代入室内温度值
yy=subs（yy,m,18）;
%确定冷却系数k
yy3=subs（yy,t,3）;
yy3=char（yy3）;
yy3=strcat（yy3,'-50=0'）;
kk=solve（yy3,k）;
double（kk）
%求10分钟后，物体的温度
yy=subs（yy,k,kk）;
yy10=subs（yy,t,10）
%求当物体的温度为30时，经过的时间
yy30=char（yy）;
yy30=strcat（yy30,'-30=0'）;
tt=solve（yy30,t）;
double（tt）
结果为：
10分钟后，物体温度为35℃左右。
当物体的温度为34℃时，经过的时间为13.8分左右。
进一步的思考
前面我们做了一些假设使问题简化，如果改变某些假设（比如说室温不是恒定不变的），我们的模型会怎么样变化。
模型中的参数若有误差，我们的结果会怎么样变化。
∙离散系统
实验目的
学会对实际问题进行数学抽象。
掌握常用组合算法和图论算法。
学会用相关软件设计程序求解。
实验原理
组合问题的基本解决方法：
1，1，从组合学基本概念、基本原理出发的所谓常规方法：
容斥原理解决计数问题
鸽笼原理解决存在性问题
递归法、母函数法解决序列变量问题
2，2，与问题所涉及的组合数学概念无关的非常规方法：
数学归纳法用天对存在性问题的分析
一一对应技术将总是转化为另一个有常规算法的问题
殊途同归法用于从不同角度讨论计数问题
数论方法把奇偶性、整除性等数论性质用于存在性问题的分析推理
图论的常用方法：
1，1，用于求解最短路问题的Dijkstra算法或标号法
2，
2，用于求最优生成树的标号法
3，3，用于求最大匹配的匈牙利算法
4，4，用于求最大流的标号法
案例：
光纤铺设问题
设有七个城市v1,v2,v3,v4,v5,v6,v7。
它们之间有通信光纤连通，其布置如下图G，且线段上数值为每条光纤的造价（单位十万元）。
问至少有几条光纤不出故障，城市间的通信才不会中断，若要造一个主干网络，这些不出故障的光纤应如何分布，才能使总的造价最低？
问题分析
这是一个求图G的生成树的问题。
因G有七个点，可知G的生成树有6条边，即至少需6条光纤不出故障，通信才不会中断。
如若造主干网，使总造价最低，则是一个求最优生成树的问题。
求最小生成树的算法（Kruskal）
（1）将G的边按权从小到大排列，不妨设为
e1,e2，…，em
（2）取T={e1}，再从e2开始依次将排好序的边加入到T中，使加入后由T导出的子图（即由T作为边集，T中的边相关联的点作为点集所确定的子图）不含圈，直至T中含有n-1条边。
注：
不构成圈的方法是标号法。
求最小生成树的Matlab程序：
function[tree,value,branch]=mintree（edge,vitex）
%这是一个求最优生成树的程序,调用格式为:
%[tree,value,branch]=mintree（edge,vitex）
%tree返回被取到的边的序号,1表示取该边,零则表示不取.
%value返回最优生成树的权值和.
%branch返回分枝数,值为1表示是树,大于1表示是森林.
%edge输入边的参数,
%第一列为各边的权;
%第二列为各边起点的序号;
%第二列为各边终点的序号.
%vitex输入图的顶点数.
%例如
%edge=[5,1,2;5,3,4,25,3,1]
%
e_num=size（edge,1）;
tree=zeros（1,e_num）;%初始化
value=0;
branch=0;
markv=zeros（1,vitex）;%给各顶点的标号
%边权排序
[eup,eorder]=sort（edge
升级会员 