第17章 有约束最优化问题
17.1
线性规划
线性规划是处理线性目标函数和线性约束的一种较为成熟的方法，目前已经广泛应用于军事、经济、工业、农业、教育、商业和社会科学等许多方面。
线性规划的求解方法主要是单纯形法（Simple Method）。该法由Dantzig 于1947年提出，以后经过多次改进。单纯形法是一种迭代算法，它从所有基本可行解的一个较小部分中通过迭代过程选出最优解。其迭代过程的一般描述如下。
① 将线性规划化为典范形式，从而可以得到一个初始基本可行解x (0) （初始顶点），将它作为迭代过程的出发点，其目标值为z (x (0) ) 。
② 寻找一个基本可行解x (1) , 使z (x (1) ) ≤z (x (0) ) 。方法是通过消去法将产生x (0) 的典范形式化为产生x (1) 的典范形式。
③ 继续寻找较好的基本可行解，使目标函数值不断改进。当某个基本可行解再也不能被其它基本可行解改进时，它就是所求的最优解。
MATLAB 优化工具箱中采用的是投影法，它是单纯形法的一种变种。 17.1.1 基本数学原理
线性规划问题的标准形式是(矩阵形式)
⎧m in z =CX --线性函数，C为常行向量⎪
AX =b --线性方程组⎨
⎪X ≥0⎩
线性规划的标准形式要求使目标函数最小化，约束条件取等式，变量x 非负。不符合
这几个条件的线性模型要首先转化成标准形式。 17.1.2 有关函数介绍
在 MA TLAB 工具箱中，可用 linprog 函数求解线性规划问题。 假设线性规划问题的数学模型为：
⎧m in ⎪⎪⎨⎪⎪⎩
f T x
Ax ≤b --线性方程组Aeq x =beq lb ≤x ≤ub
式中 f,x,b,beq,lb 和ub 为向量，A 和 Aeq 为矩阵。
linprog 函数的调用格式如下：
● x=linprog(f,A,b)——只有约束条件A*x
● x=linprog(f,A,b,Aeq,beq)——有等式约束，若没有不等式约束，此时要增加
A=[ ]、b=[ ]。
● x=linprog(f,A,b,Aeq,beq,Ib,ub)——定义设计变量x 的下界 Ib 和上界ub ，使得
x 始终在该范围内。若没有等式约束，令 Aeq=[ ]、beq=[ ]。
● x=linprog(f,A,b,Aeq,beq,Ib,ub,x0)——设置初值为x0。该选项只适用于中型问
题，默认时大型算法将忽略初值。
● x=linprog(f,A,b,Aeq,beq,Ib,ub,x0,options)——用options 指定的优化参数进行
最小化。
● [x, fvaI]=linprog(…) 返回解x 处的目标函数值fvaI 。
● [x, lambda, exitfIag]=linprog(…) 返回exitfIag 值，描述函数计算的退出条件。 ● [x, lambda, exitfIag, output]=linprog(…) 返回包含优化信息的输出变量output 。 ● [x, fval, exitfIag, output, lambda]=linprog(…) 将解x 处的拉格朗日乘子返回到
lambda 参数中。
调用格式中，lambda 参数是解x 处的拉格朗日乘子。它有以下一些属性： lambda.lower ——lambda 的下界。 lambda.upper ——lambda 的上界。
lambda.ineqlin ——lambda 的线性不等式。 lambda.eqIin ——lambda 的线性等式。 其他参数意义可参见表15-7和表15-8。
根据问题规模的不同，linprog 函数使用不同的算法：
大型优化算法——大型优化算法采用的是 LIPSOL 法。该法在进行迭代计算
之前首先要进行一系列的预处理。
中型优化算法——linprog 函数使用的是投影法，就象quadprog 函数的算法一
样。linprog 函数使用的是一种活动集方法，是线性规划中单纯形法的变种。它通过求解另一个线性规划问题来找到初始可行解。
对于大型算法，算法的第1步涉及到一些约束条件的预处理问题。有些问题可能导致linprog 函数退出，并显示不可行的消息。在本例中，exitfIag 参数将被设为负值以表示优化失败。
若 Aeq 参数中某行的所有元素都为零，但Beq 参数中对应的元素不为零，则显示以下退出消息：
Exiting due to infeasibiIity: amn aII zero row in row in the constraint matrix does not have a zero in corresponding right hand side entry.
若x 的某一个元素没在界内，则给出以下消息：
Exiting due to infeasibiIity:objectivef*x is unbounded beIow.
若Aeq 参数的某一行中只有一个非零值，则 x 中的相关值称为奇异变量。这里，x 中该成分的值可以用Aeq 和 Beq 算得。若算得的值与另一个约束条件相矛盾，则给出以下退出消息：
Exiting due to infeasibiIity:singIeton variabIes in equaIity constraits are not feasibIe.
若奇异变量可以求解但其解超出上界和下界，则给出以下退出消息：
Exiting due to infeasibiIity:singIeton variabIes in the equaIity constraints are not within
bounds.
注意：预处理步骤是累加的。例如，即使约束矩阵开始不含有元素全为零的行，其他预处理步骤也会也会引起某行元素全为零。
一旦预处理结束，将进行迭代计算，直到满足终止准则。若迭代的残值在残差在增加而不是减小，或者残差不增加也不减小，则分别给出下面两条终止消息：
One or more of the residuaIs,duaIity gap,or totaI reIative error has grown 100000 times
greater than its minimum vaIue so far: 或者
One or more of the residuaIs,duaIity gap,or totaI reIative error has staIIed:
对于中型优化问题，当解不可行时，Iinprog 函数给出下面的警告消息：
Warning: The constraints are overly stringent; there is no feasible solution. 这里，linprog 函数给出一个结果，使约束矛盾的最坏程度变到最小。 当等式约束不协调时linprog 函数给出警告消息：
Warning: The equality constraints are overly stringent; there is no feasible solution.
超出边界的解给出的警告消息是：
Warning: The soIution is unbounded and at infinity; the constraints are not restrictive
enough.
这里，linprog 函数返回 x 的值，该值满足约束条件。 另外，对于中型优化问题，显示水平参数只能使用 off 和finaI ，进行迭代输出的iter 属性不可用。
17.1.3 应用实例
【例17-1】 生产决策问题。某厂生产甲、乙两种产品，已知制成一吨产品甲需用资源A 3吨，资源B 4m 3 ；制成一吨产品乙需用资源A 2吨，资源B 6m3 ，资源C7个单位。若一吨产品甲和乙的经济价值分别为7万元和5万元，三种资源的限制量分别为90吨、200m 3 和210个单位。试决定应生产这两种产品各多少吨才能使创造的总经济价值最高?
解 令生产产品甲的数量为 x 1，生产产品乙的数量为x 2。由题意可以建立下面的模型：
⎧m ax
⎪⎪⎪⎨⎪⎪⎪⎩
z =7x 1+5x 23x 1+2x 2≤904x 1+6x 2≤200 7x 2≤210 x ≥0, x 2≥0
这是一个只有不等式约束条件的模型，该模型要求使目标函数最大化，需要按照
MATLAB 的要求进行转换，即目标函数为
m in
z =-7x 1-5x 2.
首先输入下列系数： >> f=[-7;-5];
>> A=[3 2;4 6;0 7]; >> b=[90;200;210]; >> lb=zeros(2,1);
然后调用linprog 函数； format long
>> format compact
>> [x,fval,exitflag,output,lambda]=linprog(f,A,b,[ ], [ ], lb) % 为加入lb,, 要加Aeq=[ ], beq=[ ] Optimization terminated successfully. x =
13.[1**********]554 24.[1**********]080 fval =
-2.[**************]e+002 exitflag = 1
output =
iterations: 5 cgiterations: 0
algorithm: 'lipsol' lambda =
ineqlin: [3x1 double] eqlin: [0x1 double] upper: [2x1 double] lower: [2x1 double]
由上可知，生产甲种产品14吨、乙种产品24吨可使创建的总经济价值最高。最高经济价值为218万元。exitfIag=1 表示过程正常收敛于解x 处。
【例17-2】 投资问题。某单位有一批资金用于4个工程项目的投资，各工程项目所 得到的净收益（投入资金的百分比）如表17-1所示 表17-1 工程项目收益表
由于某种原因，决定用于项目A 的投资不大于其他各项投资之和；而用于项目B 和C 的投资要大于项目D 的投资。试确定使该单位收益最大的投资分配方案。
解 用x1,x2,x3和x4分别代表用于项目 A,B,C 和D 的投资百分数，由于各项目的投资百分数之和必须等于100%，所以
x1+x2+x3+x4=1. 据题意，可以建立下面的数学模型： ⎧m ax ⎪⎪⎪
⎨
⎪⎪⎪⎩
z =0. 15x 1+0. 1x 2+0. 12x 3x 1-x 2-x 3-x 4≤0x 2+x 3-x 4≥0 x 1+x 2+x 3+x 4=1 x j ≥0, j =1, 2,..., 4
将它转换为标准形式(目标函数转化为求最小，大于等于转化为小于等于) ： ⎧m in ⎪⎪⎪
⎨
⎪⎪⎪⎩
z =-0. 15x 1-0. 1x 2-0. 12x 3x 1-x 2-x 3-x 4≤0-x 2-x 3+x 4≤0 x 1+x 2+x 3+x 4=1 x j ≥0, j =1, 2,..., 4
下面进行求解：
首先输入下列系数 >> f=[-0.15 -0.1 -0.08 -0.12]'; >> A=[1 -1 -1 -1;0 -1 -1 1]; >> b=[0;0];
>> Aeq=[1 1 1 1]; >> beq=[1];
>> lb=zeros(4,1);
然后调用Iinprog 函数： >> format short
>> [x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb) Optimization terminated successfully. x =
0.5000 0.2500 0.0000 0.2500 fval =
-0.1300 exitflag = 1 output =
iterations: 9 cgiterations: 0
algorithm: 'lipsol' lambda =
ineqlin: [2x1 double] eqlin: 0.1300
upper: [4x1 double] lower: [4x1 double]
可见，4个项目的投资百分数分别为0.50，0.25，0.00和0.25时可使该单位获得最大的收益。最大收益为13%。过程正常收敛。
【例17-3】 工件加工任务分配问题。某车间有两台机床甲和乙，可用于加工3种工件。假定这两台机床的可用台时数分别为700和800，3种工件的数量分别为300，500，和400，且已知用两台不同机床加工单位数量的不同工件所需的台时数和加工费用（如表17-2所示）。问怎样分配机床的加工任务，才能既满足加工工件的要求，又使总加工费用最低。
表17-2 机床加工情况表
解 该在甲机床上加工工件1，2和3的数量分别为x1,x2和x3，在乙机床上加工工件1，2和3的数量分别为x4,x5和x6，根据3种工件的数量限制，有 x1+x4=300 (对工件1) x2+x5=500 (对工件2) x3+x6=400 (对工件3)
再根据机床甲和乙的可用总台时限制，可以得到其他约束条件。以总加工费用最少为目标函数，组合约束条件，可以得到下面的数学模型：
⎧m in ⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩
z =13x 1+9x 2+10x 3+11x 4+12x 5+8x 6x 1+x 4=300x 2+x 5=500x 3+x 6=400
0. 4x 1+1. 1x 2+x 3≤7000. 5x 4+1. 2x 5+1. 3x 6≤800 x j ≥0, j =1, 2,..., 6
首先输入下列系数
>> f=[13 9 10 11 12 8]';
>> A=[0.4 1.1 1 0 0 0;0 0 0 0.5 1.2 1.3]; >> b=[700;800];
>> Aeq=[1 0 0 1 0 0; 0 1 0 0 1 0;0 0 1 0 0 1]; >> beq=[700;800]; >> beq=[300;500;400]; >> lb=zeros(6,1);
然后调用linprog 函数
>> [x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb) Optimization terminated successfully. x =
0.0000 500.0000 0.0000 300.0000 0.0000 400.0000 fval =
1.1000e+004 exitflag = 1 output =
iterations: 4 cgiterations: 0
algorithm: 'lipsol' lambda =
ineqlin: [2x1 double] eqlin: [3x1 double] upper: [6x1 double] lower: [6x1 double]
可见，在甲机床上加工500个工件2，在乙机床上加工300个工件1、加工400个工件3可在满足条件的情况下使总加工费最小。最小费用为11000元。收敛正常。
（以下几个例题的处理暂略，读者可以自己给出如上类似结果）
【17-4】 裁料问题。在某个建筑工程施工中需要制作10000套钢筋，每套钢筋由2.9m,2.1m 和1.5m 3种不同长度的钢筋各一根组成，它们的直径和材质不同。目前在市场
上采购到的同类钢筋的长度每根均为7.4 m ，问应购进多少根7.4m 长的钢筋才能满足工程的需要？
首先分析共有多少种不同的套裁方法，该问题的可能材料方案如表17-3所示。 表17-3 材料方案表
设以 表示按第 种裁料方案下料的原材料数量，则可得该问题的数学模型为
首先输入下列系数
然后调用Iinprog 函数：
所以最节省的情况需要9167根7.4 m 长的钢筋，其中第1种方案使用5000根，第5种方案使用1667根，第6种方案使用2500根。
【例17-5】工作人员计划安排问题。某昼夜服务的公共交通系统每天个时间段（每4小时为一个时间段）所需的值班人数如表17-4所示。这些值班人员在某一时段开始上班后要连续工作8小时（包括轮流用膳时间）。问该公交系统至少需要多少名工作人员才能满足值班的需要？
表17-4 各时段所需值班人数表
设 为第 个时段开始上班的人员数，据题意建立下面的数学模型：
需要对前面6个约束条件进行形式变换，使不等式变为非正不等式。只需要在不等式两侧取负即可。
首先输入下列系数： 然后调用Iinprog 函数：
可见，只要6个时段分别安排42人、28人、35人15人、10人、和20人就可以满足值班
的需要。共计150人。计算收敛。
【例17-6】厂址选择问题。考虑A,B,C 三地，每地都出产一定数量的原料，也消耗一定数量的产品（见表17-5）。已知制成每吨产品需3吨原料，各地之间的距离为：A~B:150km, A~C: 100km, B~C: 200km. 假定每万吨原料运输1km 的运价是5000元，每万吨产品运输1km 的运价是6000元。由于地区条件的差异，在不同地点设厂的生产费用也不同。问究竟在哪些地方设厂，规模多大，才能使总费用最小？另外，由于其他条件限制，在B 处建厂的规模（生产的产品数量）不能超过5万吨。
表17-5 A, B, C三地出产原料，消耗产品情况表 表17-5 A, B, C三地出产原料、消耗产品情况表
令 为由 地运到 地的原料数量（万吨）， 为由 地运往 地的产品数量（万吨）， （分别对应A,B,C 三地）。根据题意，可以建立问题的数学模型（其中目标函数包括原材料运输费、产品运输费和生产费）：首先输入下列系数： 然后调用Iinprog 函数：
要使总费用最小，需要B 地向A 地运送1万吨；A 、B 、C 三地的建厂规模分别为7万吨、5万吨和8万吨。最小总费用为3485元。
17.2 有约束非线性最优化问题
17.2.1 基本数学原理
在有约束最优化问题中，通常要将该问题转换为更简单的子问题，这些子问题可以求解并作为跌代过程的基础。早期的方法通常是通过构造惩罚函数等来将有约束的最优化问题转换为无约束最优化问题进行求解。现在，这些方法已经被更有效的基于K-T （Kuhn-Tucker ）方程解的方法所取代。K-T 方程是有约束最优化问题求解的必要条件。假设有所谓的Convex 规划问题，f (x ) 和G i (x ), i =1, 2,..., m 为Convex 函数，则K-T 方程对于求得全局极小点是必要的，也是充分的。 对于规划问题
⎧m in n ⎪x ∈R ⎪⎨⎪⎪⎩
f (x )
G i (x ) =0i =1,..., m e G i (x ) ≤0i =m e +1,..., m x l ≤x ≤x u
式中，x 是设计参数向量，f(x)为目标函数，返回标量值，向量函数G (x ) 返回等式约束和不等式约束在x 处的值。 它的K-T 方程可表达为
*
f (x ) +∑λ*i ⋅∇G i (x ) =0 (*)*
i =1
m
∇G i (x *) =0i =1,..., m e --梯度为零时G (?)达到极值
λ*i ≥0i =m e +1,..., m
其中第1行描述了目标函数和约束条件在解处梯度消失，需要用拉格朗日乘子
λi (i =1, 2,..., m ) 来平衡目标函数与约束梯度间大小的差异。
K-T 方程的解形成了许多非线性规划算法的基础。这些算法直接计算拉格朗日乘子。用
拟牛顿法更新过程，给K-T 方程积累二阶信息，可以保证有约束拟牛顿法的超线性收敛。这些方法称为序列二次规划法（SQP), 因为在每一次主要的迭代中都求解一次二次规划问题。 对于给定的规划问题，序列二次规划（SQP ）的主要思路是形成基于拉格朗日函数二次近似的二次规划子问题，即
L (x , λ) =f (x ) +∑λi g i (x )
i =1m
这里，通过假设约束条件为不等式约束使（*）式得到了简化，通过非线性有约束问题线性化来获得二次规划子问题。
二次规划子问题可表达为
1T
d H c d +∇f (x k ) T d d ∈R 2
∇g i (x k ) T d +g i (x k ) =0i =1,..., m e m in n
∇g i (x k ) T d +g i (x k ) ≤0i =m e +1,..., m
该子问题可以用任意一种二次规划算法求解，求得的解可以用来形成新的迭代公式
x k +1=x k +αkd k
用SQP 法求解非线性有约束问题时的迭代次数常比用无约束问题时的少，因为在搜索
区域内，SQP 方法可以获得最佳的搜索方向和步长信息。
给Rosenbrock 函数添加非线性不等式约束g(x)
22x 1+x 2-1. 5≤0
经过96次迭代得到问题的解：x=[0.9072,0.8288],初值为x=[-1.9,2],无约束问题则需要140
次迭代。图17-1是搜索路径图。
MATLAB 中 SQP 法的实现分3步，即 ①拉格朗日Hess 矩阵的更新； ②二次规划问题求解；
③一维搜索和目标函数的计算。
1．Hess 矩阵的更新
在每一次主要迭代过程中，都用BFGS 计算拉格朗日函数的Hess 矩阵的拟牛顿近似矩阵。更新公式为
H k +1=H k +
T q k q k T q k S k
-
T
H k H k T S k H k S k
式中
q k =∇f (x k +1) +∑λi ⋅∇g i (x k +1) -(∇f (x k ) +∑λi -∇g i (x k ))
i =1
i =1
n
n
上式中，λi , i =1, 2,..., m 为对拉格朗日乘子的估计。
2. 二次规划求解
SQP 法的每一次主要迭代过程都要求一次二次规划问题，形式如下：
m in q (d ) =n
d ∈R
A i d =b i A i d ≤b i
1T
d Hd +c T d 2
i =1,..., m e i =m e +1,..., m
求解过程分两步，第一步涉及可行点（若存在）的计算，第2步为可行点至解的迭代序
列。在第1步中，需要有可行点作为初值，若当前点不可行，则通过求解下列线性规划问题可以得到一个可行点：
γ∈R , x ∈R n
m in γ
i =1,..., m e
i =m e +1,..., m
A i x =b i A i x -γ≤b i
式中，A i 为矩阵A 的第i 行. 3. 一维搜索和目标函数的计算
*
二次规划子问题的解生成一个向量d k ，它形成一个新的迭代公式：
x k +1=x k +αd k
(α为步长参数)
目标函数的形式如下：
ψ(x ) =f (x ) +∑γi ⋅g i (x ) +
i =1
m e
γi ⋅m ax{0, g i (x )} ∑i =m +1
e
m
式中
γi =(γk +1) i =m ax {λi ⋅((λk ) i +λi )}i =1,..., m
i
1
2
17.2.2 有关函数介绍
利用fmincon 函数求多变量有约束非线性函数的最小值。假设多变量非线性函数的数学模型为：
m in f (x ) x
c (x ) ≤0
ceq (x ) =0 Ax ≤b
Aeq x ≤beq
lb ≤x ≤ub
式中，x ,b,beq,lb,和ub 为向量，A, 和Aeq 为矩阵，c (x) 和ceq (x ) 为函数，返回标量。f (x ), c (x ), 和ceq (x ) 可以是非线性函数。
fmincon 函数的调用格式如下：
● fmincon 求多变量有约束非线性函数的最小值，常用于有约束非线性优化问题。
● x=fmincon(fun,x0,A,b)——给定初值x0, 求解fun 函数的最小值点x 。 fun 函数的约束条
件为A*x
● x=fmincon(fun,x0,A,b,Aeq,beq)——最小化fun 函数，约束条件为Aeq*x=beq和A*x
若没有不等式存在，则设置A=[ ]，b=[ ].
● x=fmincon(fun,x0,A,b,Aeq,beq，lb ，ub ）——定义设计变量x 的下界和上界ub ，使得总
是有lb
供非线性不等式c(x)或等式ceq(x). fmincon 函数要求c(x)
● x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)——用optiions 参数指定的参数进行
最小化。
● x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2,…) ——将问题参数P1，P2，
等直接传递给函数fun 和nonlin. 若不需要这些变量，则传递空矩阵到A ，b ，Aeq ，beq ，lb ，ub, nonIcon 和options.
● [x,fval]=fmincon(…) ——返回解x 处的目标函数值。
● [x,fval,exitflag]=fmincon(…) ——返回exitfIag 参数，描述函数计算的退出条件。 ● [x,fval,exitflag,output]=fmincon(…) ——返回包含优化信息的输出参数output.
● [x,fval,exitflag,output,lambda]=fmincon(…) ——返回解x 处包含拉格朗日乘子的lambda
参数。
● [x,fval,exitflag,output,lambda,grad]=fmincon(…) ——返回解x 处fun 函数的梯度。
● [x,fval,exitflag,output,lambda,grad,hessian]=fmincon(…) ——返回解x 处fun 函数的Hess
矩阵。
各调用格式中，nonlcon 参数计算非线性不等式约束 c (x ) ≤0 和非线性等式约束ceq (x )=0。nonlcon 参数是一个包含函数名的字符串。该函数可以是M 文件、内部文件或MEX 文件。它要求输入一个向量x ，返回两个变量——解x 处的非线性不等式向量c 和非线性等式向量ceq 。
例如，若nonlcon=’mycon ’ , 则M 文件mycon.m 具有下面的形式：
function [c,ceq]=mycon(x)
c=... % 计算x 处的非线性不等式
ceq=... % 计算x 处的非线性等式
若还计算了约束的梯度，即
options=optimset('GradConstr','on')
则nonlcon 函数必须在第3个和第4个输出变量中返回c (x ) 的梯度 Gc 和ceq (x ) 的梯度Gceq 。当被调用的nonlcon 函数只需要两个输出变量（此时用户算法只需要c 和ceq 的值，而不需要Gc 和Gceq ）时，可以通过查看nargout 的值来避免计算Gc 和Gceq 的值。
function [c,ceq,Gc,Gceq]=mycon(x)
c=... % 解x 处的非线性不等式
ceq=... % 解x 处的非线性等式
if nargout > 2 % 被调用的nonlcon 函数, 要求有4个输出变量
Gc=... % 不等式的梯度
Gceq=... % 等式的梯度
end
若nonlcon 函数返回 m 元素的向量c 和长度为n 的x ，则c (x ) 的梯度Gc 是一个n ×m 的矩阵，其中 Gc (i,j ) 是c(j)对x (i ) 的偏导数。同样，若ceq 是一个 p 元素的向量，则ceq (x ) 的梯度Gceq 是一个n ×p 的矩阵，其中Gceq (i,j ) 是ceq(j ) 对x (i ) 偏导数。
其它参数意义可以参见表15-7和表15-8。
根据问题规模的不同，fmincon 函数使用不同的优化算法：
大型优化算法：默认时，若提供了函数的梯度信息，并且只有上下界存在或只有线
性等式约束存在，则fmincon 函数将选择大型算法。本法是基于内部映射牛顿法（inteior-refIective Newton method）的子空间置信域法（subspace trust-region）该法的每一次迭代都与用PCG 法求解大型线性系统得到的近似解有关。
中型优化算法：fmincon 函数使用序列二次规划法（SQP ）。本法中，在每一步迭代
中求解二次规划子问题，并用BFGS 法更新拉格朗日Hess 矩阵。
在使用该函数的过程中，还有一些需要注意的问题：
(1) 大型优化问题
①使用大型算法，必须在fun 函数中提供梯度信息(options.GradObj设置为’on ’) 。如果没有梯度信息，则给出警告消息。
fmincon 函数允许g(x)为一近似梯度，但使用真正的梯度将使优化过程更具稳健性。 ②当对矩阵的二阶导数（即Hess 矩阵）进行计算以后，用该函数求解大型问题将更有效。但不需要求得真正的Hess 矩阵，如果能提供Hess 矩阵的稀疏结构的信息（用options 参数的HessPattern 属性），则fmincon 函数可以算得Hess 矩阵的稀疏有限差分近似。
③若x0不是严格可行的，则fmincon 函数选择一个新的严格可行初始点。
④若x 的某些元素没有上界或下界，则fmincon 函数更希望对应的元素设置为Inf(对于上界) 或-Inf(对于下界），而不希望强制性地给上界赋一个很大的值，给下界赋一个很小的负值。
⑤线性约束最小化课题中也有几个问题需要注意：
Aeq 矩阵中若存在密集列或近密集列(A dense(或fairly dense)column), 将导致满
秩并使计算费时。
fmincon 函数剔除Aeq 中线性相关的行。此过程需要进行反复的因子分解，因
此，如果相关行很多的话，计算将是一件很费时的事情。
每一次迭代都要用下式进行稀疏最小二乘求解
B =Aeq T R -T
式中R T 为前提条件的乔累斯基因子。
(2) 中型优化问题
①如果用Aeq 和beq 清楚地提供等式约束，将比用lb 和ub 获得更好的数值解。
②在二次子问题中，若有等式约束并且诸因变等式(dependent equalities)被发现和剔除，则将在过程标题中显示’dependent ’(当output 参数要求使用options.Display=’iter ’）。只有在等式连续的情况下，因变等式才会被剔除。若等式系统不连续，则子问题将不可行并在过程标题中输出’infeasible ’消息。
求大型优化问题的代码中不允许上限和下限相等，即不能有lb (2)= =ub (2),否则给出下面的出错消息：
Equal upper and lower bounds not permitted in this large-scale method.
Use equality constraints and the medium-scale method instead.
若只有等式约束，仍然可以使用大型算法。当既有等式约束又有边界约束时，使用中型算法。
（3）使用fmincon 函数的一些要求
① 目标函数和约束函数都必须是连续的，否则可能会给出局部最优解。
② 当问题不可行时，fmincon 函数将试图使最大约束值最小化。
③ 目标函数和约束函数都必须是实数。
④ 对于大型优化问题，使用大型优化算法时，用户必须在fun 函数中提供梯度（options 参数的GradObj 属性必须设置为’on ’），并且只可以指定上界和下界约束，或者只有线性约束必须存在，Aeq 的行数不能多于列数。
⑤ 现在，如果在fun 函数中提供了解析梯度，则选项参数DerivativeCheck 不能与大型方法一起用，以比较解析梯度和有限差分梯度。可以通过将options 参数的MaxIter 属性设置为0来用中型方法核对导数。然后用大型方法求解问题。
17.2.3 应用实例
【例17-6】求侧面积为常数6×52㎡ 的体积最大的长方体体积。设该长方体的长、宽、高分别为x1,x2和x3.
解 据题意得到下面的数学摸型：
⎧m in z =-x 1x 2x 3 ⎨2(x x +x x +x x ) =150122331⎩
编一个M 文件, 返回x 处的函数值f 。
function f=myfun1706(x)
f=-x(1)*x(2)*x(3);
由于约束条件是非线性等式约束，所以需要编写一个约束条件M 文件:
function [c,ceq]=mycon1706(x)
c=[ ] % 书中漏了此行，害得我费了很长时间排错
ceq=x(1)*x(2)+x(2)*x(3)+x(3)*x(1)-75
下一步给定初值，并调用优化过程。
>> x0=[4;5;6];
>> lb=zeros(3,1);
>> [x,fval,exitflag,output,lambda]=fmincon(@myfun1706,x0,[],[],[],[],lb,[],@mycon1706) 计算结果为：
Warning: Large-scale (trust region) method does not currently solve this type of problem, switching to medium-scale (line search).
> In C:\MATLAB6p5\toolbox\optim\fmincon.m at line 213
Optimization terminated successfully:
First-order optimality measure less than options.TolFun and
maximum constraint violation is less than options.TolCon
Active Constraints:
1
x =
5.0000
5.0000
5.0000
fval =
-125.0000
exitflag =
1
output =
iterations: 7
funcCount: 41
stepsize: 1
algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'
firstorderopt: 4.9310e-007
cgiterations: []
lambda =
lower: [3x1 double]
upper: [3x1 double]
eqlin: [0x1 double]
eqnonlin: 2.5000
ineqlin: [0x1 double]
ineqnonlin: [0x1 double]
优化结果显示过程成功收敛，搜索方向小于两倍options.TolX, 最大违约值小于options.TolCon, 主动约束为1个。
问题的解为x (1)=x (2)=x (3)=5.0000m ，最大体积为125.0000。exitfIag=1,表示过程成功收敛于解x 处。Output 输出变量显示了收敛过程中的迭代次数、目标函数计数次数、步长、算法等信息。lambda 则包含模型信息。
【例17-7】 试设计一压缩圆柱螺旋弹簧，要求其质量最小。弹簧材料为65Mn, 最大工作载荷P =40N，最小工作载荷为0，载荷变化频率f=25Hz,弹簧寿命为104h, 弹簧钢丝直径d 的取值范围为1~4mm,中径D 的取值范围为10~30mm,工作圈数 n 不应小于4.5圈，弹簧旋绕比C 不应小于4，弹簧一端固定，一端自由，工作温度为50℃, 弹簧变形量不小于10mm 。 本题的优化目标是使弹簧质量最小，圆柱螺旋弹簧的质量可以表示为
M =γ(n +n 2) πD π
4d 2
式中：γ−弹簧材料的密度，对于钢材γ=1. 8⨯10-6kg /mm 3
N ——工作圈数；
n2——死圈数，常取n2 =1.5~2.5，现取n2=2;
D ——弹簧直径，（mm ）, 记为x(3)
d ——弹簧钢丝直径，（mm ）, 记为x(1)
n ——工作圈数，记为x(2)
将已知参数代入公式，并由
>>7.8*10^(-6)*pi*pi/4
ans =
1.9246e-005
进行整理以后得到问题的目标函数为
f (x ) =M =0. 192457⨯10-4(x (2) +2) x (1)^2x (3)
根据弹簧性能和结构上的要求，可写出问题的约束条件；
⑴强度条件
-2. 860. 86 g 1(x ) =350-163. 0x 1x 3≥0
⑵刚度条件
-43g 2(x ) =0. 4⨯10-2x 1x 2x 3-10. 0≥0
⑶稳定性条件
-43g 3(x ) =3. 7x 3-(x 2+1. 5) x 1-0. 44⨯10-2x 1x 2x 3≥0
⑷不发生共振现象，要求
-1-2 g 4(x ) =0. 356⨯106x 1x 2x 3-375≥0
⑸弹簧旋绕比的限制
-1g 5(x ) =x 3x 1-4. 0≥0
⑹对d, n, D的限制, 相当分别对x(1), x(2), x(3)的限制
1. 0≤x (1) ≤4. 0, 用应取标准值, 即1. 0, 1. 2, 1. 6, 2. 0, 2. 5, 3. 0, 3. 5, 4. 0mm 等
4. 5≤x (2) ≤50
10≤x (3) ≤30
由上可知，该压缩圆柱螺旋的优化设计是一个三维的约束优化问题，其数学模型为
min f (x ) =M =0. 192457⨯10-4(x (2) +2) x (1)^2x (3)
0. 86⎧g 1(x ) =350-163. 0x 1-2. 86x 3≤0⎪-2-43⎪g 2(x ) =0. 4⨯10x 1x 2x 3-10. 0≥0
⎪g (x ) =3. 7x -(x +1. 5) x -0. 44⨯10-2x -4x x 3≥0321123⎪3
-1-2⎪g 4(x ) =0. 356⨯106x 1x 2x 3-375≥0⎪⎪g 5(x ) =x 3/x 1-4. 0≥0 ⎪⎨g 6(x ) =x 1-1≥0
⎪g (x ) =4-x ≥01⎪7
⎪g 8(x ) =x 2-4. 5≥0⎪⎪g 9(x ) =50-x 2≥0
⎪g 10(x ) =x 3-10≥0⎪⎪g 11(x ) =30-x 3≥0⎩
取初始设计参数为 x (0) [2. 0, 5. 0, 25. 0]T .
首先编写目标函数的M 文件， 返回x 处的函数值f ：
function f=myfun1707(x)
f=0.192457*1e-4*(x(2)+2)*x(1)^2*x(3);
由于约束条件中有非线性约束，所以需要编写一个描述非线性约束条件的M 文件：
function [c,ceq]=mycon1707(x)
c(1)=350-163*x(1)^(-2.86)*x(3)^0.86;
c(2)=10-0.4*0.01*x(1)^(-4)*x(2)*x(3)^3;
c(3)=(x(2)+1.5)*x(1)+0.44*0.01*x(1)^(-4)*x(2)*x(3)^3-3.7*x(3); c(4)=375-0.356*1e6*x(1)*x(2)^(-1)*x(3)^(-2);
c(5)=4-x(3)/x(1);
ceq=[];
然后设置线性约束的系数：
>> A=[-1 0 0; 1 0 0; 0 -1 0; 0 1 0; 0 0 -1;0 0 1];
>> b=[-1 4 -4.5 50 -10 30]';
下一步给定初值，给定变量的下限约束，并调用优化过程。
>> x0=[2.0 5.0 25.0]';
>> lb=zeros(3,1);
>> [x,fval,exitflag,output,lambda]=fmincon(@myfun1707,x0,A,b,[],[],lb,[],@mycon1707) 计算结果为
Warning: Large-scale (trust region) method does not currently solve this type of problem, switching to medium-scale (line search).
> In C:\MATLAB6p5\toolbox\optim\fmincon.m at line 213
Optimization terminated successfully:
Magnitude of directional derivative in search direction
less than 2*options.TolFun and maximum constraint violation
is less than options.TolCon
Active Constraints:
6
x =
1.6564
4.5000
16.1141
fval =
0.0055
exitflag =
1
output =
iterations: 3
funcCount: 19
stepsize: 1
algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'
firstorderopt: 0.0067 % 书上这个结果为“[]”??
cgiterations: []
lambda =
lower: [3x1 double]
upper: [3x1 double]
eqlin: [0x1 double]
eqnonlin: [0x1 double]
ineqlin: [6x1 double]
ineqnonlin: [5x1 double]
所以当弹簧钢丝的直径 d 、工作圈数n 及中径D 分别取1.6564、4.5000和16.1141时弹簧质量最小，为5.5克。考虑到实际情况，各参数可分别取1.6，5.0和16.0
x=[1.6 5.0 16.0]; % 不知书上这四行从何而来？？
z=optim253(x)
z=
0.0055
故此时弹簧质量仍为5.5克。
第17章 有约束最优化问题
17.1
线性规划
线性规划是处理线性目标函数和线性约束的一种较为成熟的方法，目前已经广泛应用于军事、经济、工业、农业、教育、商业和社会科学等许多方面。
线性规划的求解方法主要是单纯形法（Simple Method）。该法由Dantzig 于1947年提出，以后经过多次改进。单纯形法是一种迭代算法，它从所有基本可行解的一个较小部分中通过迭代过程选出最优解。其迭代过程的一般描述如下。
① 将线性规划化为典范形式，从而可以得到一个初始基本可行解x (0) （初始顶点），将它作为迭代过程的出发点，其目标值为z (x (0) ) 。
② 寻找一个基本可行解x (1) , 使z (x (1) ) ≤z (x (0) ) 。方法是通过消去法将产生x (0) 的典范形式化为产生x (1) 的典范形式。
③ 继续寻找较好的基本可行解，使目标函数值不断改进。当某个基本可行解再也不能被其它基本可行解改进时，它就是所求的最优解。
MATLAB 优化工具箱中采用的是投影法，它是单纯形法的一种变种。 17.1.1 基本数学原理
线性规划问题的标准形式是(矩阵形式)
⎧m in z =CX --线性函数，C为常行向量⎪
AX =b --线性方程组⎨
⎪X ≥0⎩
线性规划的标准形式要求使目标函数最小化，约束条件取等式，变量x 非负。不符合
这几个条件的线性模型要首先转化成标准形式。 17.1.2 有关函数介绍
在 MA TLAB 工具箱中，可用 linprog 函数求解线性规划问题。 假设线性规划问题的数学模型为：
⎧m in ⎪⎪⎨⎪⎪⎩
f T x
Ax ≤b --线性方程组Aeq x =beq lb ≤x ≤ub
式中 f,x,b,beq,lb 和ub 为向量，A 和 Aeq 为矩阵。
linprog 函数的调用格式如下：
● x=linprog(f,A,b)——只有约束条件A*x
● x=linprog(f,A,b,Aeq,beq)——有等式约束，若没有不等式约束，此时要增加
A=[ ]、b=[ ]。
● x=linprog(f,A,b,Aeq,beq,Ib,ub)——定义设计变量x 的下界 Ib 和上界ub ，使得
x 始终在该范围内。若没有等式约束，令 Aeq=[ ]、beq=[ ]。
● x=linprog(f,A,b,Aeq,beq,Ib,ub,x0)——设置初值为x0。该选项只适用于中型问
题，默认时大型算法将忽略初值。
● x=linprog(f,A,b,Aeq,beq,Ib,ub,x0,options)——用options 指定的优化参数进行
最小化。
● [x, fvaI]=linprog(…) 返回解x 处的目标函数值fvaI 。
● [x, lambda, exitfIag]=linprog(…) 返回exitfIag 值，描述函数计算的退出条件。 ● [x, lambda, exitfIag, output]=linprog(…) 返回包含优化信息的输出变量output 。 ● [x, fval, exitfIag, output, lambda]=linprog(…) 将解x 处的拉格朗日乘子返回到
lambda 参数中。
调用格式中，lambda 参数是解x 处的拉格朗日乘子。它有以下一些属性： lambda.lower ——lambda 的下界。 lambda.upper ——lambda 的上界。
lambda.ineqlin ——lambda 的线性不等式。 lambda.eqIin ——lambda 的线性等式。 其他参数意义可参见表15-7和表15-8。
根据问题规模的不同，linprog 函数使用不同的算法：
大型优化算法——大型优化算法采用的是 LIPSOL 法。该法在进行迭代计算
之前首先要进行一系列的预处理。
中型优化算法——linprog 函数使用的是投影法，就象quadprog 函数的算法一
样。linprog 函数使用的是一种活动集方法，是线性规划中单纯形法的变种。它通过求解另一个线性规划问题来找到初始可行解。
对于大型算法，算法的第1步涉及到一些约束条件的预处理问题。有些问题可能导致linprog 函数退出，并显示不可行的消息。在本例中，exitfIag 参数将被设为负值以表示优化失败。
若 Aeq 参数中某行的所有元素都为零，但Beq 参数中对应的元素不为零，则显示以下退出消息：
Exiting due to infeasibiIity: amn aII zero row in row in the constraint matrix does not have a zero in corresponding right hand side entry.
若x 的某一个元素没在界内，则给出以下消息：
Exiting due to infeasibiIity:objectivef*x is unbounded beIow.
若Aeq 参数的某一行中只有一个非零值，则 x 中的相关值称为奇异变量。这里，x 中该成分的值可以用Aeq 和 Beq 算得。若算得的值与另一个约束条件相矛盾，则给出以下退出消息：
Exiting due to infeasibiIity:singIeton variabIes in equaIity constraits are not feasibIe.
若奇异变量可以求解但其解超出上界和下界，则给出以下退出消息：
Exiting due to infeasibiIity:singIeton variabIes in the equaIity constraints are not within
bounds.
注意：预处理步骤是累加的。例如，即使约束矩阵开始不含有元素全为零的行，其他预处理步骤也会也会引起某行元素全为零。
一旦预处理结束，将进行迭代计算，直到满足终止准则。若迭代的残值在残差在增加而不是减小，或者残差不增加也不减小，则分别给出下面两条终止消息：
One or more of the residuaIs,duaIity gap,or totaI reIative error has grown 100000 times
greater than its minimum vaIue so far: 或者
One or more of the residuaIs,duaIity gap,or totaI reIative error has staIIed:
对于中型优化问题，当解不可行时，Iinprog 函数给出下面的警告消息：
Warning: The constraints are overly stringent; there is no feasible solution. 这里，linprog 函数给出一个结果，使约束矛盾的最坏程度变到最小。 当等式约束不协调时linprog 函数给出警告消息：
Warning: The equality constraints are overly stringent; there is no feasible solution.
超出边界的解给出的警告消息是：
Warning: The soIution is unbounded and at infinity; the constraints are not restrictive
enough.
这里，linprog 函数返回 x 的值，该值满足约束条件。 另外，对于中型优化问题，显示水平参数只能使用 off 和finaI ，进行迭代输出的iter 属性不可用。
17.1.3 应用实例
【例17-1】 生产决策问题。某厂生产甲、乙两种产品，已知制成一吨产品甲需用资源A 3吨，资源B 4m 3 ；制成一吨产品乙需用资源A 2吨，资源B 6m3 ，资源C7个单位。若一吨产品甲和乙的经济价值分别为7万元和5万元，三种资源的限制量分别为90吨、200m 3 和210个单位。试决定应生产这两种产品各多少吨才能使创造的总经济价值最高?
解 令生产产品甲的数量为 x 1，生产产品乙的数量为x 2。由题意可以建立下面的模型：
⎧m ax
⎪⎪⎪⎨⎪⎪⎪⎩
z =7x 1+5x 23x 1+2x 2≤904x 1+6x 2≤200 7x 2≤210 x ≥0, x 2≥0
这是一个只有不等式约束条件的模型，该模型要求使目标函数最大化，需要按照
MATLAB 的要求进行转换，即目标函数为
m in
z =-7x 1-5x 2.
首先输入下列系数： >> f=[-7;-5];
>> A=[3 2;4 6;0 7]; >> b=[90;200;210]; >> lb=zeros(2,1);
然后调用linprog 函数； format long
>> format compact
>> [x,fval,exitflag,output,lambda]=linprog(f,A,b,[ ], [ ], lb) % 为加入lb,, 要加Aeq=[ ], beq=[ ] Optimization terminated successfully. x =
13.[1**********]554 24.[1**********]080 fval =
-2.[**************]e+002 exitflag = 1
output =
iterations: 5 cgiterations: 0
algorithm: 'lipsol' lambda =
ineqlin: [3x1 double] eqlin: [0x1 double] upper: [2x1 double] lower: [2x1 double]
由上可知，生产甲种产品14吨、乙种产品24吨可使创建的总经济价值最高。最高经济价值为218万元。exitfIag=1 表示过程正常收敛于解x 处。
【例17-2】 投资问题。某单位有一批资金用于4个工程项目的投资，各工程项目所 得到的净收益（投入资金的百分比）如表17-1所示 表17-1 工程项目收益表
由于某种原因，决定用于项目A 的投资不大于其他各项投资之和；而用于项目B 和C 的投资要大于项目D 的投资。试确定使该单位收益最大的投资分配方案。
解 用x1,x2,x3和x4分别代表用于项目 A,B,C 和D 的投资百分数，由于各项目的投资百分数之和必须等于100%，所以
x1+x2+x3+x4=1. 据题意，可以建立下面的数学模型： ⎧m ax ⎪⎪⎪
⎨
⎪⎪⎪⎩
z =0. 15x 1+0. 1x 2+0. 12x 3x 1-x 2-x 3-x 4≤0x 2+x 3-x 4≥0 x 1+x 2+x 3+x 4=1 x j ≥0, j =1, 2,..., 4
将它转换为标准形式(目标函数转化为求最小，大于等于转化为小于等于) ： ⎧m in ⎪⎪⎪
⎨
⎪⎪⎪⎩
z =-0. 15x 1-0. 1x 2-0. 12x 3x 1-x 2-x 3-x 4≤0-x 2-x 3+x 4≤0 x 1+x 2+x 3+x 4=1 x j ≥0, j =1, 2,..., 4
下面进行求解：
首先输入下列系数 >> f=[-0.15 -0.1 -0.08 -0.12]'; >> A=[1 -1 -1 -1;0 -1 -1 1]; >> b=[0;0];
>> Aeq=[1 1 1 1]; >> beq=[1];
>> lb=zeros(4,1);
然后调用Iinprog 函数： >> format short
>> [x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb) Optimization terminated successfully. x =
0.5000 0.2500 0.0000 0.2500 fval =
-0.1300 exitflag = 1 output =
iterations: 9 cgiterations: 0
algorithm: 'lipsol' lambda =
ineqlin: [2x1 double] eqlin: 0.1300
upper: [4x1 double] lower: [4x1 double]
可见，4个项目的投资百分数分别为0.50，0.25，0.00和0.25时可使该单位获得最大的收益。最大收益为13%。过程正常收敛。
【例17-3】 工件加工任务分配问题。某车间有两台机床甲和乙，可用于加工3种工件。假定这两台机床的可用台时数分别为700和800，3种工件的数量分别为300，500，和400，且已知用两台不同机床加工单位数量的不同工件所需的台时数和加工费用（如表17-2所示）。问怎样分配机床的加工任务，才能既满足加工工件的要求，又使总加工费用最低。
表17-2 机床加工情况表
解 该在甲机床上加工工件1，2和3的数量分别为x1,x2和x3，在乙机床上加工工件1，2和3的数量分别为x4,x5和x6，根据3种工件的数量限制，有 x1+x4=300 (对工件1) x2+x5=500 (对工件2) x3+x6=400 (对工件3)
再根据机床甲和乙的可用总台时限制，可以得到其他约束条件。以总加工费用最少为目标函数，组合约束条件，可以得到下面的数学模型：
⎧m in ⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩
z =13x 1+9x 2+10x 3+11x 4+12x 5+8x 6x 1+x 4=300x 2+x 5=500x 3+x 6=400
0. 4x 1+1. 1x 2+x 3≤7000. 5x 4+1. 2x 5+1. 3x 6≤800 x j ≥0, j =1, 2,..., 6
首先输入下列系数
>> f=[13 9 10 11 12 8]';
>> A=[0.4 1.1 1 0 0 0;0 0 0 0.5 1.2 1.3]; >> b=[700;800];
>> Aeq=[1 0 0 1 0 0; 0 1 0 0 1 0;0 0 1 0 0 1]; >> beq=[700;800]; >> beq=[300;500;400]; >> lb=zeros(6,1);
然后调用linprog 函数
>> [x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb) Optimization terminated successfully. x =
0.0000 500.0000 0.0000 300.0000 0.0000 400.0000 fval =
1.1000e+004 exitflag = 1 output =
iterations: 4 cgiterations: 0
algorithm: 'lipsol' lambda =
ineqlin: [2x1 double] eqlin: [3x1 double] upper: [6x1 double] lower: [6x1 double]
可见，在甲机床上加工500个工件2，在乙机床上加工300个工件1、加工400个工件3可在满足条件的情况下使总加工费最小。最小费用为11000元。收敛正常。
（以下几个例题的处理暂略，读者可以自己给出如上类似结果）
【17-4】 裁料问题。在某个建筑工程施工中需要制作10000套钢筋，每套钢筋由2.9m,2.1m 和1.5m 3种不同长度的钢筋各一根组成，它们的直径和材质不同。目前在市场
上采购到的同类钢筋的长度每根均为7.4 m ，问应购进多少根7.4m 长的钢筋才能满足工程的需要？
首先分析共有多少种不同的套裁方法，该问题的可能材料方案如表17-3所示。 表17-3 材料方案表
设以 表示按第 种裁料方案下料的原材料数量，则可得该问题的数学模型为
首先输入下列系数
然后调用Iinprog 函数：
所以最节省的情况需要9167根7.4 m 长的钢筋，其中第1种方案使用5000根，第5种方案使用1667根，第6种方案使用2500根。
【例17-5】工作人员计划安排问题。某昼夜服务的公共交通系统每天个时间段（每4小时为一个时间段）所需的值班人数如表17-4所示。这些值班人员在某一时段开始上班后要连续工作8小时（包括轮流用膳时间）。问该公交系统至少需要多少名工作人员才能满足值班的需要？
表17-4 各时段所需值班人数表
设 为第 个时段开始上班的人员数，据题意建立下面的数学模型：
需要对前面6个约束条件进行形式变换，使不等式变为非正不等式。只需要在不等式两侧取负即可。
首先输入下列系数： 然后调用Iinprog 函数：
可见，只要6个时段分别安排42人、28人、35人15人、10人、和20人就可以满足值班
的需要。共计150人。计算收敛。
【例17-6】厂址选择问题。考虑A,B,C 三地，每地都出产一定数量的原料，也消耗一定数量的产品（见表17-5）。已知制成每吨产品需3吨原料，各地之间的距离为：A~B:150km, A~C: 100km, B~C: 200km. 假定每万吨原料运输1km 的运价是5000元，每万吨产品运输1km 的运价是6000元。由于地区条件的差异，在不同地点设厂的生产费用也不同。问究竟在哪些地方设厂，规模多大，才能使总费用最小？另外，由于其他条件限制，在B 处建厂的规模（生产的产品数量）不能超过5万吨。
表17-5 A, B, C三地出产原料，消耗产品情况表 表17-5 A, B, C三地出产原料、消耗产品情况表
令 为由 地运到 地的原料数量（万吨）， 为由 地运往 地的产品数量（万吨）， （分别对应A,B,C 三地）。根据题意，可以建立问题的数学模型（其中目标函数包括原材料运输费、产品运输费和生产费）：首先输入下列系数： 然后调用Iinprog 函数：
要使总费用最小，需要B 地向A 地运送1万吨；A 、B 、C 三地的建厂规模分别为7万吨、5万吨和8万吨。最小总费用为3485元。
17.2 有约束非线性最优化问题
17.2.1 基本数学原理
在有约束最优化问题中，通常要将该问题转换为更简单的子问题，这些子问题可以求解并作为跌代过程的基础。早期的方法通常是通过构造惩罚函数等来将有约束的最优化问题转换为无约束最优化问题进行求解。现在，这些方法已经被更有效的基于K-T （Kuhn-Tucker ）方程解的方法所取代。K-T 方程是有约束最优化问题求解的必要条件。假设有所谓的Convex 规划问题，f (x ) 和G i (x ), i =1, 2,..., m 为Convex 函数，则K-T 方程对于求得全局极小点是必要的，也是充分的。 对于规划问题
⎧m in n ⎪x ∈R ⎪⎨⎪⎪⎩
f (x )
G i (x ) =0i =1,..., m e G i (x ) ≤0i =m e +1,..., m x l ≤x ≤x u
式中，x 是设计参数向量，f(x)为目标函数，返回标量值，向量函数G (x ) 返回等式约束和不等式约束在x 处的值。 它的K-T 方程可表达为
*
f (x ) +∑λ*i ⋅∇G i (x ) =0 (*)*
i =1
m
∇G i (x *) =0i =1,..., m e --梯度为零时G (?)达到极值
λ*i ≥0i =m e +1,..., m
其中第1行描述了目标函数和约束条件在解处梯度消失，需要用拉格朗日乘子
λi (i =1, 2,..., m ) 来平衡目标函数与约束梯度间大小的差异。
K-T 方程的解形成了许多非线性规划算法的基础。这些算法直接计算拉格朗日乘子。用
拟牛顿法更新过程，给K-T 方程积累二阶信息，可以保证有约束拟牛顿法的超线性收敛。这些方法称为序列二次规划法（SQP), 因为在每一次主要的迭代中都求解一次二次规划问题。 对于给定的规划问题，序列二次规划（SQP ）的主要思路是形成基于拉格朗日函数二次近似的二次规划子问题，即
L (x , λ) =f (x ) +∑λi g i (x )
i =1m
这里，通过假设约束条件为不等式约束使（*）式得到了简化，通过非线性有约束问题线性化来获得二次规划子问题。
二次规划子问题可表达为
1T
d H c d +∇f (x k ) T d d ∈R 2
∇g i (x k ) T d +g i (x k ) =0i =1,..., m e m in n
∇g i (x k ) T d +g i (x k ) ≤0i =m e +1,..., m
该子问题可以用任意一种二次规划算法求解，求得的解可以用来形成新的迭代公式
x k +1=x k +αkd k
用SQP 法求解非线性有约束问题时的迭代次数常比用无约束问题时的少，因为在搜索
区域内，SQP 方法可以获得最佳的搜索方向和步长信息。
给Rosenbrock 函数添加非线性不等式约束g(x)
22x 1+x 2-1. 5≤0
经过96次迭代得到问题的解：x=[0.9072,0.8288],初值为x=[-1.9,2],无约束问题则需要140
次迭代。图17-1是搜索路径图。
MATLAB 中 SQP 法的实现分3步，即 ①拉格朗日Hess 矩阵的更新； ②二次规划问题求解；
③一维搜索和目标函数的计算。
1．Hess 矩阵的更新
在每一次主要迭代过程中，都用BFGS 计算拉格朗日函数的Hess 矩阵的拟牛顿近似矩阵。更新公式为
H k +1=H k +
T q k q k T q k S k
-
T
H k H k T S k H k S k
式中
q k =∇f (x k +1) +∑λi ⋅∇g i (x k +1) -(∇f (x k ) +∑λi -∇g i (x k ))
i =1
i =1
n
n
上式中，λi , i =1, 2,..., m 为对拉格朗日乘子的估计。
2. 二次规划求解
SQP 法的每一次主要迭代过程都要求一次二次规划问题，形式如下：
m in q (d ) =n
d ∈R
A i d =b i A i d ≤b i
1T
d Hd +c T d 2
i =1,..., m e i =m e +1,..., m
求解过程分两步，第一步涉及可行点（若存在）的计算，第2步为可行点至解的迭代序
列。在第1步中，需要有可行点作为初值，若当前点不可行，则通过求解下列线性规划问题可以得到一个可行点：
γ∈R , x ∈R n
m in γ
i =1,..., m e
i =m e +1,..., m
A i x =b i A i x -γ≤b i
式中，A i 为矩阵A 的第i 行. 3. 一维搜索和目标函数的计算
*
二次规划子问题的解生成一个向量d k ，它形成一个新的迭代公式：
x k +1=x k +αd k
(α为步长参数)
目标函数的形式如下：
ψ(x ) =f (x ) +∑γi ⋅g i (x ) +
i =1
m e
γi ⋅m ax{0, g i (x )} ∑i =m +1
e
m
式中
γi =(γk +1) i =m ax {λi ⋅((λk ) i +λi )}i =1,..., m
i
1
2
17.2.2 有关函数介绍
利用fmincon 函数求多变量有约束非线性函数的最小值。假设多变量非线性函数的数学模型为：
m in f (x ) x
c (x ) ≤0
ceq (x ) =0 Ax ≤b
Aeq x ≤beq
lb ≤x ≤ub
式中，x ,b,beq,lb,和ub 为向量，A, 和Aeq 为矩阵，c (x) 和ceq (x ) 为函数，返回标量。f (x ), c (x ), 和ceq (x ) 可以是非线性函数。
fmincon 函数的调用格式如下：
● fmincon 求多变量有约束非线性函数的最小值，常用于有约束非线性优化问题。
● x=fmincon(fun,x0,A,b)——给定初值x0, 求解fun 函数的最小值点x 。 fun 函数的约束条
件为A*x
● x=fmincon(fun,x0,A,b,Aeq,beq)——最小化fun 函数，约束条件为Aeq*x=beq和A*x
若没有不等式存在，则设置A=[ ]，b=[ ].
● x=fmincon(fun,x0,A,b,Aeq,beq，lb ，ub ）——定义设计变量x 的下界和上界ub ，使得总
是有lb
供非线性不等式c(x)或等式ceq(x). fmincon 函数要求c(x)
● x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)——用optiions 参数指定的参数进行
最小化。
● x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2,…) ——将问题参数P1，P2，
等直接传递给函数fun 和nonlin. 若不需要这些变量，则传递空矩阵到A ，b ，Aeq ，beq ，lb ，ub, nonIcon 和options.
● [x,fval]=fmincon(…) ——返回解x 处的目标函数值。
● [x,fval,exitflag]=fmincon(…) ——返回exitfIag 参数，描述函数计算的退出条件。 ● [x,fval,exitflag,output]=fmincon(…) ——返回包含优化信息的输出参数output.
● [x,fval,exitflag,output,lambda]=fmincon(…) ——返回解x 处包含拉格朗日乘子的lambda
参数。
● [x,fval,exitflag,output,lambda,grad]=fmincon(…) ——返回解x 处fun 函数的梯度。
● [x,fval,exitflag,output,lambda,grad,hessian]=fmincon(…) ——返回解x 处fun 函数的Hess
矩阵。
各调用格式中，nonlcon 参数计算非线性不等式约束 c (x ) ≤0 和非线性等式约束ceq (x )=0。nonlcon 参数是一个包含函数名的字符串。该函数可以是M 文件、内部文件或MEX 文件。它要求输入一个向量x ，返回两个变量——解x 处的非线性不等式向量c 和非线性等式向量ceq 。
例如，若nonlcon=’mycon ’ , 则M 文件mycon.m 具有下面的形式：
function [c,ceq]=mycon(x)
c=... % 计算x 处的非线性不等式
ceq=... % 计算x 处的非线性等式
若还计算了约束的梯度，即
options=optimset('GradConstr','on')
则nonlcon 函数必须在第3个和第4个输出变量中返回c (x ) 的梯度 Gc 和ceq (x ) 的梯度Gceq 。当被调用的nonlcon 函数只需要两个输出变量（此时用户算法只需要c 和ceq 的值，而不需要Gc 和Gceq ）时，可以通过查看nargout 的值来避免计算Gc 和Gceq 的值。
function [c,ceq,Gc,Gceq]=mycon(x)
c=... % 解x 处的非线性不等式
ceq=... % 解x 处的非线性等式
if nargout > 2 % 被调用的nonlcon 函数, 要求有4个输出变量
Gc=... % 不等式的梯度
Gceq=... % 等式的梯度
end
若nonlcon 函数返回 m 元素的向量c 和长度为n 的x ，则c (x ) 的梯度Gc 是一个n ×m 的矩阵，其中 Gc (i,j ) 是c(j)对x (i ) 的偏导数。同样，若ceq 是一个 p 元素的向量，则ceq (x ) 的梯度Gceq 是一个n ×p 的矩阵，其中Gceq (i,j ) 是ceq(j ) 对x (i ) 偏导数。
其它参数意义可以参见表15-7和表15-8。
根据问题规模的不同，fmincon 函数使用不同的优化算法：
大型优化算法：默认时，若提供了函数的梯度信息，并且只有上下界存在或只有线
性等式约束存在，则fmincon 函数将选择大型算法。本法是基于内部映射牛顿法（inteior-refIective Newton method）的子空间置信域法（subspace trust-region）该法的每一次迭代都与用PCG 法求解大型线性系统得到的近似解有关。
中型优化算法：fmincon 函数使用序列二次规划法（SQP ）。本法中，在每一步迭代
中求解二次规划子问题，并用BFGS 法更新拉格朗日Hess 矩阵。
在使用该函数的过程中，还有一些需要注意的问题：
(1) 大型优化问题
①使用大型算法，必须在fun 函数中提供梯度信息(options.GradObj设置为’on ’) 。如果没有梯度信息，则给出警告消息。
fmincon 函数允许g(x)为一近似梯度，但使用真正的梯度将使优化过程更具稳健性。 ②当对矩阵的二阶导数（即Hess 矩阵）进行计算以后，用该函数求解大型问题将更有效。但不需要求得真正的Hess 矩阵，如果能提供Hess 矩阵的稀疏结构的信息（用options 参数的HessPattern 属性），则fmincon 函数可以算得Hess 矩阵的稀疏有限差分近似。
③若x0不是严格可行的，则fmincon 函数选择一个新的严格可行初始点。
④若x 的某些元素没有上界或下界，则fmincon 函数更希望对应的元素设置为Inf(对于上界) 或-Inf(对于下界），而不希望强制性地给上界赋一个很大的值，给下界赋一个很小的负值。
⑤线性约束最小化课题中也有几个问题需要注意：
Aeq 矩阵中若存在密集列或近密集列(A dense(或fairly dense)column), 将导致满
秩并使计算费时。
fmincon 函数剔除Aeq 中线性相关的行。此过程需要进行反复的因子分解，因
此，如果相关行很多的话，计算将是一件很费时的事情。
每一次迭代都要用下式进行稀疏最小二乘求解
B =Aeq T R -T
式中R T 为前提条件的乔累斯基因子。
(2) 中型优化问题
①如果用Aeq 和beq 清楚地提供等式约束，将比用lb 和ub 获得更好的数值解。
②在二次子问题中，若有等式约束并且诸因变等式(dependent equalities)被发现和剔除，则将在过程标题中显示’dependent ’(当output 参数要求使用options.Display=’iter ’）。只有在等式连续的情况下，因变等式才会被剔除。若等式系统不连续，则子问题将不可行并在过程标题中输出’infeasible ’消息。
求大型优化问题的代码中不允许上限和下限相等，即不能有lb (2)= =ub (2),否则给出下面的出错消息：
Equal upper and lower bounds not permitted in this large-scale method.
Use equality constraints and the medium-scale method instead.
若只有等式约束，仍然可以使用大型算法。当既有等式约束又有边界约束时，使用中型算法。
（3）使用fmincon 函数的一些要求
① 目标函数和约束函数都必须是连续的，否则可能会给出局部最优解。
② 当问题不可行时，fmincon 函数将试图使最大约束值最小化。
③ 目标函数和约束函数都必须是实数。
④ 对于大型优化问题，使用大型优化算法时，用户必须在fun 函数中提供梯度（options 参数的GradObj 属性必须设置为’on ’），并且只可以指定上界和下界约束，或者只有线性约束必须存在，Aeq 的行数不能多于列数。
⑤ 现在，如果在fun 函数中提供了解析梯度，则选项参数DerivativeCheck 不能与大型方法一起用，以比较解析梯度和有限差分梯度。可以通过将options 参数的MaxIter 属性设置为0来用中型方法核对导数。然后用大型方法求解问题。
17.2.3 应用实例
【例17-6】求侧面积为常数6×52㎡ 的体积最大的长方体体积。设该长方体的长、宽、高分别为x1,x2和x3.
解 据题意得到下面的数学摸型：
⎧m in z =-x 1x 2x 3 ⎨2(x x +x x +x x ) =150122331⎩
编一个M 文件, 返回x 处的函数值f 。
function f=myfun1706(x)
f=-x(1)*x(2)*x(3);
由于约束条件是非线性等式约束，所以需要编写一个约束条件M 文件:
function [c,ceq]=mycon1706(x)
c=[ ] % 书中漏了此行，害得我费了很长时间排错
ceq=x(1)*x(2)+x(2)*x(3)+x(3)*x(1)-75
下一步给定初值，并调用优化过程。
>> x0=[4;5;6];
>> lb=zeros(3,1);
>> [x,fval,exitflag,output,lambda]=fmincon(@myfun1706,x0,[],[],[],[],lb,[],@mycon1706) 计算结果为：
Warning: Large-scale (trust region) method does not currently solve this type of problem, switching to medium-scale (line search).
> In C:\MATLAB6p5\toolbox\optim\fmincon.m at line 213
Optimization terminated successfully:
First-order optimality measure less than options.TolFun and
maximum constraint violation is less than options.TolCon
Active Constraints:
1
x =
5.0000
5.0000
5.0000
fval =
-125.0000
exitflag =
1
output =
iterations: 7
funcCount: 41
stepsize: 1
algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'
firstorderopt: 4.9310e-007
cgiterations: []
lambda =
lower: [3x1 double]
upper: [3x1 double]
eqlin: [0x1 double]
eqnonlin: 2.5000
ineqlin: [0x1 double]
ineqnonlin: [0x1 double]
优化结果显示过程成功收敛，搜索方向小于两倍options.TolX, 最大违约值小于options.TolCon, 主动约束为1个。
问题的解为x (1)=x (2)=x (3)=5.0000m ，最大体积为125.0000。exitfIag=1,表示过程成功收敛于解x 处。Output 输出变量显示了收敛过程中的迭代次数、目标函数计数次数、步长、算法等信息。lambda 则包含模型信息。
【例17-7】 试设计一压缩圆柱螺旋弹簧，要求其质量最小。弹簧材料为65Mn, 最大工作载荷P =40N，最小工作载荷为0，载荷变化频率f=25Hz,弹簧寿命为104h, 弹簧钢丝直径d 的取值范围为1~4mm,中径D 的取值范围为10~30mm,工作圈数 n 不应小于4.5圈，弹簧旋绕比C 不应小于4，弹簧一端固定，一端自由，工作温度为50℃, 弹簧变形量不小于10mm 。 本题的优化目标是使弹簧质量最小，圆柱螺旋弹簧的质量可以表示为
M =γ(n +n 2) πD π
4d 2
式中：γ−弹簧材料的密度，对于钢材γ=1. 8⨯10-6kg /mm 3
N ——工作圈数；
n2——死圈数，常取n2 =1.5~2.5，现取n2=2;
D ——弹簧直径，（mm ）, 记为x(3)
d ——弹簧钢丝直径，（mm ）, 记为x(1)
n ——工作圈数，记为x(2)
将已知参数代入公式，并由
>>7.8*10^(-6)*pi*pi/4
ans =
1.9246e-005
进行整理以后得到问题的目标函数为
f (x ) =M =0. 192457⨯10-4(x (2) +2) x (1)^2x (3)
根据弹簧性能和结构上的要求，可写出问题的约束条件；
⑴强度条件
-2. 860. 86 g 1(x ) =350-163. 0x 1x 3≥0
⑵刚度条件
-43g 2(x ) =0. 4⨯10-2x 1x 2x 3-10. 0≥0
⑶稳定性条件
-43g 3(x ) =3. 7x 3-(x 2+1. 5) x 1-0. 44⨯10-2x 1x 2x 3≥0
⑷不发生共振现象，要求
-1-2 g 4(x ) =0. 356⨯106x 1x 2x 3-375≥0
⑸弹簧旋绕比的限制
-1g 5(x ) =x 3x 1-4. 0≥0
⑹对d, n, D的限制, 相当分别对x(1), x(2), x(3)的限制
1. 0≤x (1) ≤4. 0, 用应取标准值, 即1. 0, 1. 2, 1. 6, 2. 0, 2. 5, 3. 0, 3. 5, 4. 0mm 等
4. 5≤x (2) ≤50
10≤x (3) ≤30
由上可知，该压缩圆柱螺旋的优化设计是一个三维的约束优化问题，其数学模型为
min f (x ) =M =0. 192457⨯10-4(x (2) +2) x (1)^2x (3)
0. 86⎧g 1(x ) =350-163. 0x 1-2. 86x 3≤0⎪-2-43⎪g 2(x ) =0. 4⨯10x 1x 2x 3-10. 0≥0
⎪g (x ) =3. 7x -(x +1. 5) x -0. 44⨯10-2x -4x x 3≥0321123⎪3
-1-2⎪g 4(x ) =0. 356⨯106x 1x 2x 3-375≥0⎪⎪g 5(x ) =x 3/x 1-4. 0≥0 ⎪⎨g 6(x ) =x 1-1≥0
⎪g (x ) =4-x ≥01⎪7
⎪g 8(x ) =x 2-4. 5≥0⎪⎪g 9(x ) =50-x 2≥0
⎪g 10(x ) =x 3-10≥0⎪⎪g 11(x ) =30-x 3≥0⎩
取初始设计参数为 x (0) [2. 0, 5. 0, 25. 0]T .
首先编写目标函数的M 文件， 返回x 处的函数值f ：
function f=myfun1707(x)
f=0.192457*1e-4*(x(2)+2)*x(1)^2*x(3);
由于约束条件中有非线性约束，所以需要编写一个描述非线性约束条件的M 文件：
function [c,ceq]=mycon1707(x)
c(1)=350-163*x(1)^(-2.86)*x(3)^0.86;
c(2)=10-0.4*0.01*x(1)^(-4)*x(2)*x(3)^3;
c(3)=(x(2)+1.5)*x(1)+0.44*0.01*x(1)^(-4)*x(2)*x(3)^3-3.7*x(3); c(4)=375-0.356*1e6*x(1)*x(2)^(-1)*x(3)^(-2);
c(5)=4-x(3)/x(1);
ceq=[];
然后设置线性约束的系数：
>> A=[-1 0 0; 1 0 0; 0 -1 0; 0 1 0; 0 0 -1;0 0 1];
>> b=[-1 4 -4.5 50 -10 30]';
下一步给定初值，给定变量的下限约束，并调用优化过程。
>> x0=[2.0 5.0 25.0]';
>> lb=zeros(3,1);
>> [x,fval,exitflag,output,lambda]=fmincon(@myfun1707,x0,A,b,[],[],lb,[],@mycon1707) 计算结果为
Warning: Large-scale (trust region) method does not currently solve this type of problem, switching to medium-scale (line search).
> In C:\MATLAB6p5\toolbox\optim\fmincon.m at line 213
Optimization terminated successfully:
Magnitude of directional derivative in search direction
less than 2*options.TolFun and maximum constraint violation
is less than options.TolCon
Active Constraints:
6
x =
1.6564
4.5000
16.1141
fval =
0.0055
exitflag =
1
output =
iterations: 3
funcCount: 19
stepsize: 1
algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'
firstorderopt: 0.0067 % 书上这个结果为“[]”??
cgiterations: []
lambda =
lower: [3x1 double]
upper: [3x1 double]
eqlin: [0x1 double]
eqnonlin: [0x1 double]
ineqlin: [6x1 double]
ineqnonlin: [5x1 double]
所以当弹簧钢丝的直径 d 、工作圈数n 及中径D 分别取1.6564、4.5000和16.1141时弹簧质量最小，为5.5克。考虑到实际情况，各参数可分别取1.6，5.0和16.0
x=[1.6 5.0 16.0]; % 不知书上这四行从何而来？？
z=optim253(x)
z=
0.0055
故此时弹簧质量仍为5.5克。