This slide is based on the book:
赵静、但琦 主编《数学建模与数学实验》（第4版）
在优化模型中, 如果目标函数或约束条件中至少有一个是非线性函数, 则称此最优化问题为
非线性规划问题的数学模型一般可以表示为 \[ \begin{aligned} &\min f(X)\\ &\text{s.t.}\ \begin{cases} g_i(X)\geqslant 0,\quad i=1,2,\ldots,m\\ h_j(X)=0,\quad i=1,2,\ldots,\ell\\ \end{cases} \end{aligned} \]
其中 $X=(x_1,x_2,\ldots,x_n)^T\in\mathbb{E}^n$, $f,g_i,h_j$ 是定义在 $\mathbb{E}^n$ 上的实值函数, 简记为 \[ \begin{aligned} f:\ \mathbb{E}^n\rightarrow\mathbb{E}^1\\ g_i:\ \mathbb{E}^n\rightarrow\mathbb{E}^1\\ h_j:\ \mathbb{E}^n\rightarrow\mathbb{E}^1\\ \end{aligned} \]
如果采用向量表示法, 则可以写为 \[ \begin{aligned} &\min f(X)\\ &\text{s.t.}\ \begin{cases} g(X)\geqslant 0,\\ h(X)=0. \end{cases} \end{aligned} \]
其中 $g(X)=(g_1(X),g_2(X),\ldots,g_m(X))^T$, $h(X)=(h_1(X),h_2(X),\ldots,h_\ell(X))^T$. 即 $g,h$ 分别定义在 $\mathbb{E}^n$ 上而取值于 $\mathbb{E}^m$ 和 $\mathbb{E}^{\ell}$ 的向量值函数, 简记为 \[ \begin{aligned} g:\ \mathbb{E}^n\rightarrow\mathbb{E}^m,\\ h:\ \mathbb{E}^n\rightarrow\mathbb{E}^{\ell},\\ \end{aligned} \]
没有约束的优化问题称为
某市燃气公司计划要建一个煤气供应站, 该站要向城市中有固定位置的 $m$ 个用户供货, 对于选定的坐标系而言, 已知第 $i$ 个用户的位置(坐标)为 $P_i=(a_i,b_i)$, $i=1,2,\ldots,m$, 如果只考虑直线距离, 问如何确定这个煤气站的位置, 才能使总的运输距离最短?
设煤气站的位置为 $P=(x,y)$, 则第 $i$ 个用户到该站的直线距离为 \[ d(P,P_i)=\sqrt{(x-a_i)^2+(y-b_i)^2}, \] 故 $m$ 个用户到该站的总距离为 \[ f(P)=\sum_{i=1}^{m}d(P,P_i)=\sum_{i=1}^{m}\sqrt{(x-a_i)^2+(y-b_i)^2}. \]
选址问题归结为求点 $P\in\Omega$, 使得 $f(P)$ 达到最小值.
某大学希望为毕业生安排工作位置. 为了简单起见, 假设每个毕业生接受政府部门(Government Department)、工业界(Industry)或科学院(Academy of Sciences)中的一个位置.
令 $N_j$ 表示第 $j$ 年毕业的人数 $(j=1,2,\ldots,n)$, $G_j$, $I_j$ 和 $A_j$ 分别表示第 $j$ 年进入政府部门、工业界或科学院的人数. 即有 \[G_j+I_j+A_j=N_j.\]
考虑的一种方法是假设给出每年学生进入各行业的比例. 如果进入政府部门、工业界或科学院的比例分别为 $\lambda_1,\lambda_2,\lambda_3$, 则在第 $j$ 年可估计出进入各行业的人数为 \[ \begin{aligned} \hat{G}_j&=\lambda_1 N_j,\\ \hat{I}_j&=\lambda_2 N_j,\\ \hat{A}_j&=\lambda_3 N_j.\\ \end{aligned} \]
为了有根据地衡量这个模型的可靠性, 必须了解进入这三个行业的实际人数 $G_j,I_j,A_j$ 与预估数字 $\hat{G}_j$, $\hat{I}_j$, $\hat{A}_j$ 之间的差别.
按最小二乘法估计, 为使 \[ \sum_{j=1}^{n}\Bigl[(G_j-\hat{G}_j)^2+(I_j-\hat{I}_j)^2+(A_j-\hat{A}_j)^2\Bigr] \] 最小, 同时满足所有毕业生为这些行业所雇佣的约束条件. 根据进入各行业的比例, 模型可表示为 \[ \begin{aligned} &\min\sum_{j=1}^{n}\Bigl[(G_j-\lambda_1 N_j)^2+(I_j-\lambda_2 N_j)^2+(A_j-\lambda_3 N_j)^2\Bigr]\\ &\text{s.t.}\ \begin{cases} \lambda_1+\lambda_2+\lambda_3=1,\\ \lambda_i\geqslant 0,\quad i=1,2,3. \end{cases} \end{aligned} \]
假设某公司在下一个计划期内可用于投资的总资本为 $b$ 万元, 可供选择的投资项目共有 $n$ 个, 分别记为 $1,2,\ldots,n$. 已知对第 $j$ 个项目的投资总额为 $a_j$ 万元, 而收益总额为 $c_j$ 万元. 问如何进行投资, 才能使利润率(即单位投资可获得的收益, i.e., 收益/投资)最高?
设投资决策变量为 \[ x_j=\begin{cases} 1, & \text{若对第j个项目投资},\\ 0, & \text{若对第j个项目不投资} \end{cases} \]
则该问题可归结为求变量 $x_j(j=1,2,\ldots,n)$, 使 \[ \max f(x_1,x_2,\ldots,x_n)=\dfrac{\sum_{j=1}^{n}c_j x_j}{\sum_{j=1}^{n}a_j x_j} \] \[ \text{s.t.}\ \begin{cases} \sum_{j=1}^{n}a_j x_j\leqslant b,\\ x_j=0 \text{或}\ 1\quad(j=1,2,\ldots,n) \end{cases} \]
\[ \min F(X) \] \[ \text{s.t.}\ \begin{cases} AX\leqslant b,\\%linear A_{eq}\cdot X=b_{eq},\\%linear G(X)\leqslant 0,\\%non-linear C_{eq}(X)=0,\\%non-linear VLB\leqslant X\leqslant VUB\\ \end{cases} \]
这里 $X\in\mathbb{R}^n$, 及 $n$ 维实变元向量. $G(X)$ 与 $C_{eq}(X)$ 均为非线性函数组成的向量. 其他变量的含义与线性规划、二次规划中相同.
分为以下三个基本步骤
function f=fun(X); f=F(X);
function [G,Ceq]=nonlcon(X); G=... Ceq=...
x=fmincon('fun',X0,A,b)
x=fmincon('fun',X0,A,b,Aeq,beq)
x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB)
x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB,'nonlcon')
x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB,'nonlcon',options)
[x,fval]=fmincon(...)
[x,fval,exitflag]=fmincon(...)
[x,fval,exitflag,output]=fmincon(...)
>> help fmincon
fmincon finds a constrained minimum of a function of several variables.
fmincon attempts to solve problems of the form:
min F(X) subject to: A*X <= B, Aeq*X = Beq (linear constraints)
X C(X) <= 0, Ceq(X) = 0 (nonlinear constraints)
LB <= X <= UB (bounds)
fmincon implements four different algorithms: interior point, SQP,
active set, and trust region reflective. Choose one via the option
Algorithm: for instance, to choose SQP, set OPTIONS =
optimoptions('fmincon','Algorithm','sqp'), and then pass OPTIONS to
fmincon.
X = fmincon(FUN,X0,A,B) starts at X0 and finds a minimum X to the
function FUN, subject to the linear inequalities A*X <= B. FUN accepts
input X and returns a scalar function value F evaluated at X. X0 may be
a scalar, vector, or matrix.
X = fmincon(FUN,X0,A,B,Aeq,Beq) minimizes FUN subject to the linear
equalities Aeq*X = Beq as well as A*X <= B. (Set A=[] and B=[] if no
inequalities exist.)
X = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB) defines a set of lower and upper
bounds on the design variables, X, so that a solution is found in
the range LB <= X <= UB. Use empty matrices for LB and UB
if no bounds exist. Set LB(i) = -Inf if X(i) is unbounded below;
set UB(i) = Inf if X(i) is unbounded above.
X = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON) subjects the minimization
to the constraints defined in NONLCON. The function NONLCON accepts X
and returns the vectors C and Ceq, representing the nonlinear
inequalities and equalities respectively. fmincon minimizes FUN such
that C(X) <= 0 and Ceq(X) = 0. (Set LB = [] and/or UB = [] if no bounds
exist.)
X = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS) minimizes with
the default optimization parameters replaced by values in OPTIONS, an
argument created with the OPTIMOPTIONS function. See OPTIMOPTIONS for
details. For a list of options accepted by fmincon refer to the
documentation.
X = fmincon(PROBLEM) finds the minimum for PROBLEM. PROBLEM is a
structure with the function FUN in PROBLEM.objective, the start point
in PROBLEM.x0, the linear inequality constraints in PROBLEM.Aineq
and PROBLEM.bineq, the linear equality constraints in PROBLEM.Aeq and
PROBLEM.beq, the lower bounds in PROBLEM.lb, the upper bounds in
PROBLEM.ub, the nonlinear constraint function in PROBLEM.nonlcon, the
options structure in PROBLEM.options, and solver name 'fmincon' in
PROBLEM.solver. Use this syntax to solve at the command line a problem
exported from OPTIMTOOL.
[X,FVAL] = fmincon(FUN,X0,...) returns the value of the objective
function FUN at the solution X.
[X,FVAL,EXITFLAG] = fmincon(FUN,X0,...) returns an EXITFLAG that
describes the exit condition. Possible values of EXITFLAG and the
corresponding exit conditions are listed below. See the documentation
for a complete description.
All algorithms:
1 First order optimality conditions satisfied.
0 Too many function evaluations or iterations.
-1 Stopped by output/plot function.
-2 No feasible point found.
Trust-region-reflective, interior-point, and sqp:
2 Change in X too small.
Trust-region-reflective:
3 Change in objective function too small.
Active-set only:
4 Computed search direction too small.
5 Predicted change in objective function too small.
Interior-point and sqp:
-3 Problem seems unbounded.
[X,FVAL,EXITFLAG,OUTPUT] = fmincon(FUN,X0,...) returns a structure
OUTPUT with information such as total number of iterations, and final
objective function value. See the documentation for a complete list.
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA] = fmincon(FUN,X0,...) returns the
Lagrange multipliers at the solution X: LAMBDA.lower for LB,
LAMBDA.upper for UB, LAMBDA.ineqlin is for the linear inequalities,
LAMBDA.eqlin is for the linear equalities, LAMBDA.ineqnonlin is for the
nonlinear inequalities, and LAMBDA.eqnonlin is for the nonlinear
equalities.
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD] = fmincon(FUN,X0,...) returns the
value of the gradient of FUN at the solution X.
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = fmincon(FUN,X0,...)
returns the value of the exact or approximate Hessian of the Lagrangian
at X.
Examples
FUN can be specified using @:
X = fmincon(@humps,...)
In this case, F = humps(X) returns the scalar function value F of
the HUMPS function evaluated at X.
FUN can also be an anonymous function:
X = fmincon(@(x) 3*sin(x(1))+exp(x(2)),[1;1],[],[],[],[],[0 0])
returns X = [0;0].
If FUN or NONLCON are parameterized, you can use anonymous functions to
capture the problem-dependent parameters. Suppose you want to minimize
the objective given in the function myfun, subject to the nonlinear
constraint mycon, where these two functions are parameterized by their
second argument a1 and a2, respectively. Here myfun and mycon are
MATLAB file functions such as
function f = myfun(x,a1)
f = x(1)^2 + a1*x(2)^2;
function [c,ceq] = mycon(x,a2)
c = a2/x(1) - x(2);
ceq = [];
To optimize for specific values of a1 and a2, first assign the values
to these two parameters. Then create two one-argument anonymous
functions that capture the values of a1 and a2, and call myfun and
mycon with two arguments. Finally, pass these anonymous functions to
fmincon:
a1 = 2; a2 = 1.5; % define parameters first
options = optimoptions('fmincon','Algorithm','interior-point'); % run interior-point algorithm
x = fmincon(@(x) myfun(x,a1),[1;2],[],[],[],[],[],[],@(x) mycon(x,a2),options)
See also optimoptions, optimtool, fminunc, fminbnd, fminsearch, @, function_handle.
\[ \min f=-x_1-2x_2+\frac{1}{2}x_1^2+\frac{1}{2}x_2^2 \] \[ \text{s.t.}\ \begin{cases} 2x_1+3x_2\leqslant 6,\\ x_1+4x_2\leqslant 5,\\ x_1,x_2\geqslant 0. \end{cases} \]
写成标准型 \[ \min f=-x_1-2x_2+\frac{1}{2}x_1^2+\frac{1}{2}x_2^2 \] \[ \text{s.t.}\ \begin{cases} \begin{pmatrix} 2 & 3\\ 1 & 4\\ \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ \end{pmatrix}\leqslant \begin{pmatrix} 6\\ 5\\ \end{pmatrix}\\ \begin{pmatrix} 0\\ 0\\ \end{pmatrix}\leqslant \begin{pmatrix} x_1\\ x_2\\ \end{pmatrix} \end{cases} \]
先建立 M 文件 sec4_2_exmp1.m
% sec4_2_exmp1 % page 55 function f=sec4_2_exmp1(x) f = -x(1)-2*x(2)+(1/2)*x(1)^2+(1/2)*x(2)^2;
然后建立主程序 youh2.m (youh 是优化的拼音缩写)
% youh2.m
% together with sec4_2_exmp1.m
% page 55
x0=[1;1];
A=[2 3;1 4];
b=[6;5];
Aeq=[];
beq=[];
VLB=[0;0];
VUB=[];
[x,fval]=fmincon('sec4_2_exmp1',x0,A,b,Aeq,beq,VLB,VUB)
>> sec4_2_exmp1_youh Local minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the default value of the optimality tolerance, and constraints are satisfied to within the default value of the constraint tolerance. <stopping criteria details>
\[ \min f(X)=e^{x_1}(4x_1^2+2x_2^2+4x_1x_2+2x_2+1) \] \[ \text{s.t.}\ \begin{cases} x_1+x_2=0,\\ 1.5+x_1x_2-x_1-x_2\leqslant 0,\\ -x_1x_2-10\leqslant 0. \end{cases} \]
建立 sec4_2_exmp2.m 文件, 定义目标函数:
% sec4_2_exmp2
% page 55
% $\min f(x)=e^{x_1}(4x_1^2+2x_2^2+4x_1x_2+2x_2+1)$.
function f=sec4_2_exmp2(x)
f = exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
再建立 sec4_2_exmp2_mycon.m 文件, 定义非线性约束:
function[g,ceq]=sec4_2_exmp2_mycon(x) g=[1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10]; ceq=[];
最后建立主程序
% sec4_2_exmp2_youh.m
% together with sec4_2_exmp2.m
% page 55
x0=[-1;1];
A=[];
b=[];
Aeq=[1 1];
beq=[0];
VLB=[];
VUB=[];
[x,fval]=fmincon('sec4_2_exmp2',x0,A,b,Aeq,beq,VLB,VUB,'sec4_2_exmp2_mycon')
>> sec4_2_exmp2_youh
Local minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
x =
-3.1623
3.1623
fval =
1.1566
注: 与书上的结果不一致.
\[ \min f(X)=-2x_1-x_2 \] \[ \text{s.t.}\ \begin{cases} g_1(X)=25-x_1^2-x_2^2\geqslant 0\Rightarrow x_1^2+x_2^2-25\leqslant 0,\\ g_2(X)=7-x_1^2+x_2^2\geqslant 0 \Rightarrow x_1^2-x_2^2-7\leqslant 0,\\ 0\leqslant x_1\leqslant 5,\quad 0\leqslant x_2\leqslant 10. \end{cases} \]
建立 sec4_2_exmp3.m 文件, 定义目标函数:
% sec4_2_exmp3 % page 56 % $\min f(x)=-2x_1-x_2$. function f=sec4_2_exmp3(x) f = -2*x(1)-x(2);
再建立 sec4_2_exmp3_mycon.m 文件, 定义非线性约束:
% sec4_2_exmp3_mycon.m
% page 56
% $\min f(X)=-2x_1-x_2$.
% -----------------------
% Constrains
% --------------
% s.t.
% \begin{cases}
% g_1(X)&=25-x_1^2-x_2^2\geqslant 0,\\
% g_2(X)&=7-x_1^2+x_2^2\geqslant 0,\\
% 0\leqslant x_1\leqslant 5,\ 0\leqslant x_2\leqslant 10.
% \end{cases}
function[g,ceq]=sec4_2_exmp3_mycon(x)
g=[x(1)^2+x(2)^2-25;x(1)^2-x(2)^2-7];
ceq=[];
最后建立主程序 sec4_2_exmp3_youh.m
% sec4_2_exmp3_youh.m
% together with sec4_2_exmp3.m
% page 56
x0=[3;2.5];
A=[];
b=[];
Aeq=[];
beq=[];
VLB=[0;0];
VUB=[5;10];
[x,fval,exitflag,output]=fmincon('sec4_2_exmp3',x0,A,b,Aeq,beq,VLB,VUB,'sec4_2_exmp3_mycon')
>> sec4_2_exmp3_youh
Local minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
x =
4.0000
3.0000
fval =
-11.0000
exitflag =
1
output =
包含以下字段的 struct:
iterations: 8
funcCount: 27
constrviolation: 0
stepsize: 4.8567e-06
algorithm: 'interior-point'
firstorderopt: 4.1645e-08
cgiterations: 0
message: 'Local minimum found that satisfies the constraints.…'
某公司有6个建筑工地要开工, 每个工地的位置(用平面坐标 $(a,b)$ 表示, 距离单位: km)及水泥日用量 $d$ (单位: t)由下表给出.
| 1 | 2 | 3 | 4 | 5 | 6 | |
|---|---|---|---|---|---|---|
| a | 1.25 | 8.75 | 0.5 | 5.75 | 3 | 7.25 |
| b | 1.25 | 0.75 | 4.75 | 5 | 6.5 | 7.75 |
| d | 3 | 5 | 4 | 7 | 6 | 11 |
记工地的位置为 $(a_i,b_i)$, 水泥日用量为 $d_i$, $i=1,2,3,4,5,6$. 料场位置为 $(x_j,y_j)$, 日储量为 $e_j$, $j=1,2$.
从料场 $j$ 向工地 $i$ 运送水泥量为 $X_{ij}$ 吨.
因此, 目标函数为 \[ \min f=\sum_{j=1}^{2}\sum_{i=1}^{6}X_{ij}\cdot\sqrt{(x_j-a_i)^2+(y_j-b_i)^2} \]
约束条件为 \[ \begin{cases} \sum_{j=1}^{2}X_{ij}=d_i,\quad i=1,2,\ldots,6.\\ \sum_{i=1}^{6}X_{ij}\leqslant e_j,\quad j=1,2. \end{cases} \]
使用两个临时料场 $A(5,1)$, $B(2,7)$ ($j=1$ 对应料场 $A$, $j=2$ 对应料场 $B$). 求从料场 $j$ 向工地 $i$ 的运送量 $X_{ij}$, 在各工地用量必须满足, 以及各料场运送量不超过日储量的条件下, 使得总的吨千米数最小. 此时是一个线性规划问题.
线性规划模型为: \[ \min f=\sum_{j=1}^{2}\sum_{i=1}^{6}aa{i,j}X_{ij} \] \[ \text{s.t.}\ \begin{cases} \sum_{j=1}^{2}X_{ij}=d_i,\quad i=1,2,\ldots,6.\\ \sum_{i=1}^{6}X_{ij}\leqslant e_j,\quad j=1,2. \end{cases} \]
其中 $aa(i,j)=\sqrt{(x_j-a_i)^2+(y_j-b_i)^2}$, $i=1,2,\ldots,6$, $j=1,2$.
设 \[ \vec{x}= \begin{pmatrix} x_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\\ x_7\\ x_8\\ x_9\\ x_{10}\\ x_{11}\\ x_{12}\\ \end{pmatrix}= \begin{pmatrix} X_{11}\\ X_{21}\\ X_{31}\\ X_{41}\\ X_{51}\\ X_{61}\\ X_{12}\\ X_{22}\\ X_{32}\\ X_{42}\\ X_{52}\\ X_{62}\\ \end{pmatrix} \]
% sec4_2_exmp4
% page 57-58
% 使用临时料场
clear
a=[1.25 8.75 0.5 5.75 3 7.25];%
b=[1.25 0.75 4.75 5 6.5 7.75];
d=[3 5 4 7 6 11];
x=[5 2];
y=[1 7];
e=[20 20];
aa=zeros(6,2);%预分配内存给aa矩阵.
for i=1:6
for j=1:2
aa(i,j)=sqrt((x(j)-a(i))^2+(y(j)-b(i))^2);
end
end
% aa 矩阵是6x2型矩阵, 将其第一列和第二列并为一列，然后再转置成为行向量,
% 作为线性规划问题中的目标函数的系数向量.
CC=[aa(:,1);aa(:,2)];
A=[1 1 1 1 1 1 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 1 1];
b=[20;20];
Aeq=[1 0 0 0 0 0 1 0 0 0 0 0
0 1 0 0 0 0 0 1 0 0 0 0
0 0 1 0 0 0 0 0 1 0 0 0
0 0 0 1 0 0 0 0 0 1 0 0
0 0 0 0 1 0 0 0 0 0 1 0
0 0 0 0 0 1 0 0 0 0 0 1 ];
%beq=[d(1);d(2);d(3);d(4);d(5);d(6);];%直接写成 beq=d'
beq=d';
VLB=[0 0 0 0 0 0 0 0 0 0 0 0]';
VUB=[];
x0=[1 2 3 0 1 0 0 1 0 1 0 1]';
[x, fval]=linprog(CC,A,b,Aeq,beq,VLB,VUB,x0)
>> sec4_2_exmp4_1
警告: Your current settings will run a
different algorithm ('dual-simplex') in a
future release.
> In linprog (line 204)
In sec4_2_exmp4_1 (line 36)
The interior-point-legacy algorithm uses a built-in starting point;
ignoring supplied X0.
Optimization terminated.
x =
3.0000
5.0000
0.0000
7.0000
0.0000
1.0000
0.0000
0.0000
4.0000
0.0000
6.0000
10.0000
fval =
136.2275
由料场 $A,B$ 向6个工地运输水泥的方案为:
| 1 | 2 | 3 | 4 | 5 | 6 | |
|---|---|---|---|---|---|---|
| 料场 A | 3 | 5 | 0 | 7 | 0 | 1 |
| 料场 B | 0 | 0 | 4 | 0 | 6 | 10 |
改建两个新料场, 要同时确定料场的位置 $(x_j,y_j)$, $(j=1,2)$ 和运送量 $X_{ij}$, 在同样条件下使总的吨千米数最小. 此时是非线性规划问题.
非线性规划模型为: \[ \min f=\sum_{j=1}^{2}\sum_{i=1}^{6}X_{ij}\cdot\sqrt{(x_j-a_i)^2+(y_j-b_i)^2} \] \[ \text{s.t.}\ \begin{cases} \sum_{j=1}^{2}X_{ij}=d_i,\quad i=1,2,\ldots,6.\\ \sum_{i=1}^{6}X_{ij}\leqslant e_j,\quad j=1,2. \end{cases} \]
设 \[ \vec{x}= \begin{pmatrix} x_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\\ x_7\\ x_8\\ x_9\\ x_{10}\\ x_{11}\\ x_{12}\\ x_{13}\\ x_{14}\\ x_{15}\\ x_{16}\\ \end{pmatrix}= \begin{pmatrix} X_{11}\\ X_{21}\\ X_{31}\\ X_{41}\\ X_{51}\\ X_{61}\\ X_{12}\\ X_{22}\\ X_{32}\\ X_{42}\\ X_{52}\\ X_{62}\\ x_1\\ y_1\\ x_2\\ y_2\\ \end{pmatrix} \]
% sec4_2_exmp4_2_liaochang.m
% page 58-60
% 改建料场
function f=sec4_2_exmp4_2_liaochang(x)
a=[1.25 8.75 0.5 5.75 3 7.25];
b=[1.25 0.75 4.75 5 6.5 7.75];
d=[3 5 4 7 6 11];
e=[20 20];
f1=0;
s=zeros(1,12);
for i=1:6
s(i)=sqrt((x(13)-a(i))^2+(x(14)-b(i))^2);
f1=s(i)*x(i)+f1;
end
f2=0;
for i=7:12
s(i)=sqrt((x(15)-a(i-6))^2+(x(16)-b(i-6))^2);
f2=s(i)*x(i)+f2;
end
f=f1+f2;
end
取初值为线性规划的计算结果以及临时料场的坐标所构成的向量.
x0=[3 5 0 7 0 1 0 0 4 0 6 10 | 5 1 2 7]';
% sec4_2_exmp4_2_gying.m
% page 58-60
% 改建料场
clear
%初始料场的位置 A:(x1,y1)=(5,1) B:(x2,y2)=(2,7)
x0=[3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7]';
A=[1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0];
b=[20;20];
%Aeq 是 6x16 型矩阵.
Aeq=[1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 ];
beq=[3 5 4 7 6 11]';
VLB=[zeros(12,1); -inf; -inf; -inf; -inf];
VUB=[];
[x, fval, exitflag]=fmincon('sec4_2_exmp4_2_liaochang',x0,A,b,Aeq,beq,VLB,VUB)
%-----------
%运行结果
% x =
% 2.9410
% 4.8404
% 3.8779
% 6.9431
% 1.3034
% 0.0221
% 0.0590
% 0.1596
% 0.1221
% 0.0569
% 4.6966
% 10.9779
% 5.7298
% 4.9757
% 7.2500
% 7.7500
% fval =
% 90.4921
% exitflag =
% 2
这个结果和书上的不一致, 如果用书上计算出的 $\vec{x}$, 我们调用这里定义的 sec4_2_exmp4_2_liaochang() 函数, 结果是 89.8835
>> x2=[3 5 4 7 1 0 0 0 0 0 5 11 5.6959 4.9285 7.2500 7.7500]';
>> x2
x2 =
3.0000
5.0000
4.0000
7.0000
1.0000
0
0
0
0
0
5.0000
11.0000
5.6959
4.9285
7.2500
7.7500
>> size(x2)
ans =
16 1
>> sec4_2_exmp4_2_liaochang(x2)
ans =
89.8835
某炼油厂将 4 种不同含硫量的液体原料（分别记为甲、乙、丙、丁）混合生产两种产品（分别记为 $A,B$）.
按照生产工艺的要求, 原料甲、乙、丁必须首先倒入混合池中混合, 混合后的液体再分别与原料丙混合生产 $A,B$.
已知原料甲、乙、丙、丁的含硫量分别是 3%, 1%, 2%, 1%. 进货价格分别是 6千元/吨, 16千元/吨, 10千元/吨, 15千元/吨.
产品 $A,B$ 的含硫量分别不超过 2.5%, 1.5%. 售价分别是 9千元/吨, 15千元/吨.
根据市场信息, 原料甲、乙的供应没有限制, 原料丙、丁的供应量最多为 250吨, 100吨.
产品 $A,B$ 的市场需求量分别为 300吨, 500吨.
问应该如何安排生产?
假设液体原料在混合生产过程中没有质量损失.
根据题意, 炼油厂需要合理安排生产, 目的是使利润最大化, 故优化目标是总利润最大.
\[ \max\ \Bigl[(9-6x_1-16x_2-15x_4)\cdot y_1+(15-6x_1-16x_2-15x_4)\cdot y_2+(9-10)z_1+(15-10)z_2\Bigr] \]
原料最大供应限制 \[ \begin{aligned} x_4(y_1+y_2)\leqslant 100,\\ z_1+z_2\leqslant 250. \end{aligned} \]
产品最大需求量限制 \[ \begin{aligned} y_1+z_1\leqslant 300,\\ y_2+z_2\leqslant 500. \end{aligned} \]
产品 $A$ 的含硫量限制 \[ \frac{(3x_1+x_2+x_4)\cdot y_1+2z_1}{y_1+z_1}\leqslant 2.5 \]
产品 $B$ 的含硫量限制 \[ \frac{(3x_1+x_2+x_4)\cdot y_2+2z_2}{y_2+z_2}\leqslant 1.5 \]
比例之和为 1 \[ x_1+x_2+x_4=1 \]
非负约束 \[ \begin{aligned} x_1,x_2,x_4\in[0,1]\\ y_1,z_1\in[0,300]\\ y_2,z_2\in[0,500]\\ \end{aligned} \]
综上所述, 炼油厂生产安排问题的数学模型为:
\[ \max\ \Bigl[(9-6x_1-16x_2-15x_4)\cdot y_1+(15-6x_1-16x_2-15x_4)\cdot y_2-z_1+5z_2\Bigr] \]
\[ \text{s.t.}\ \begin{cases} z_1+z_2\leqslant 250,\\ y_1+z_1\leqslant 300,\\ y_2+z_2\leqslant 500,\\ x_1+x_2+x_4=1,\\ x_4(y_1+y_2)\leqslant 100,\\ (3x_1+x_2+x_4-2.5)\cdot y_1-0.5z_1\leqslant 0,\\ (3x_1+x_2+x_4-1.5)\cdot y_2+0.5z_2\leqslant 0,\\ 0\leqslant x_1,x_2,x_4\leqslant 1,\\ 0\leqslant y_1,z_1\leqslant 300,\\ 0\leqslant y_2,z_2\leqslant 500.\\ \end{cases} \]
这是一个非线性规划模型.
为便于求解, 将模型化为标准型
\[ \min\ -\Bigl[(9-6x_1-16x_2-15x_4)\cdot y_1+(15-6x_1-16x_2-15x_4)\cdot y_2-z_1+5z_2\Bigr] \]
\[ \text{s.t.}\ \begin{cases} z_1+z_2\leqslant 250,\\ y_1+z_1\leqslant 300,\\ y_2+z_2\leqslant 500,\\ x_1+x_2+x_4=1,\\ x_4(y_1+y_2)\leqslant 100,\\ (3x_1+x_2+x_4-2.5)\cdot y_1-0.5z_1\leqslant 0,\\ (3x_1+x_2+x_4-1.5)\cdot y_2+0.5z_2\leqslant 0,\\ 0\leqslant x_1,x_2,x_4\leqslant 1,\\ 0\leqslant y_1,z_1\leqslant 300,\\ 0\leqslant y_2,z_2\leqslant 500.\\ \end{cases} \]
为便于编程, 令 \[ x(1)=x_1,\ x(2)=x_2,\ x(3)=x_4,\ x(4)=y_1,\ x(5)=y_2,\ x(6)=z_1,\ x(7)=z_2. \]