非线性规划模型

非线性规划模型

概述

​ 如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。一般来说，解非线性规划要比线性规划问题困难得多。而且也不像线性规划有单纯形法这一通用方法，非线性规划目前还没有适用于各种问题的一般算法，各个方法都有自己特定的适用范围。

数学模型

例1 投资决策问题

​ 某企业有$\text{n}$个项目可供选择投资，并且至少要对其中一个项目投资。已知该企业拥有总资金$\text{A}$元，投资于第$\text{i}(\text{i}=1,2,\cdots ,\text{n})$个项目需要花资金$a_{\text{i}}$元，并预计可收益${\text{b}}_{\text{i}}$元。式选择最佳投资方案。

​ 解：设投资决策变量为：
$$x_{\text{i}}=\{\begin{array}{ll}1& 决定会投资第\text{i}个项目\\ 0& 决定不投资第\text{i}个项目\end{array}$$
​ 则投资总额为$\sum _{\text{i}=1}^{\text{n}}{a}_{\text{i}}{x}_{\text{i}}$，投资总收益为$\sum _{\text{i}=1}^{\text{n}}{\text{b}}_{\text{i}}{x}_{\text{i}}$。

​ 因为该公司至少要对一个项目投资，并且总的投资金额不能超过总资金$\text{A}$，固有限制条件：
$$0<\sum _{\text{i}=1}^{\text{n}}{a}_{\text{i}}{x}_{\text{i}}\le \text{A}$$
​ 另外由于$x_{\text{i}}$只能取0或1，所以还有：
$$x_{\text{i}}(1-x_{\text{i}})=0\text{i}=1,2,\cdots ,\text{n}$$
​ 最佳方案应是投资额最小而总收益最大的方案，所以这个最佳投资决策问题归结为总资金以及决策变量的限制条件下，极大化总收益和总投资之比。因此，其数学模型为：
$$\max z=\frac{\sum _{\text{i}=1}^{\text{n}}{\text{b}}_{\text{i}}{x}_{\text{i}}}{\sum _{\text{i}=1}^{\text{n}}{a}_{\text{i}}{x}_{\text{i}}}$$

$$\text{s.t.}\{\begin{array}{l}0<\sum _{\text{i}=1}^{\text{n}}{a}_{\text{i}}{x}_{\text{i}}\le \text{A}\\ x_{\text{i}}(1-x_{\text{i}})=0\text{i}=1,2,\cdots ,\text{n}\end{array}$$

非线性数学规划的数学模型一般形式：
$$\min f(x)$$

$$\text{s.t.}\{\begin{array}{ll}{h}_{\text{j}}(x)\le 0& \text{j}=1,2,\cdots ,\text{q}\\ g_{\text{i}}(x)=0& \text{i}=1,2,\cdots ,\text{p}\end{array}$$

​ 在一组等式或不等式的约束下，求一个函数的最大值（或最小值）问题，其中至少有一个非线性函数，这类问题称之为非线性规划问题。

$\text{MATLAB}$中的数学模型

一般形式

$$\min f(x)$$

$$\text{s.t.}\{\begin{array}{l}\begin{array}{r}\text{A}x\le \text{b}\\ \text{c}(x)\le 0\\ \text{ceq}(x)=0\\ \text{lb}\le x\le \text{ub}\end{array}\end{array}$$

​ 其中$f(x)$是标量函数，$\text{A},\text{b},\text{Aeq},\text{beq},\text{lb},\text{ub}$是相应维数的矩阵和向量，$\text{c}(x),\text{ceq}(x)$是非线性向量函数。

求解命令

[x,fval]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nolcon,options);

​ x的返回值是决策向量$\mathbf{X}$的取值，fval返回的是目标函数的取值，其中fun是用$\text{M}$文件定义的函数；x0是$\mathbf{X}$的初始值（初始值不同，则最优解可能不同）；A,b,Aeq,beq定义了线性约束$\text{A}\mathbf{x}\le \text{b},\text{A}\mathbf{eq}\cdot \mathbf{x}=\text{b}\mathbf{eq}$，如果没有线性约束，则A=[],b=[],Aeq=[],beq=[]；lb和ub是变量$x$的下界和上界，如果上界和下界没有约束，即$x$无下界也无上界，则lb=[],ub=[]，也可以写成lb的各分量都是-inf，ub的各分量都为inf；nonlcon是用$\text{M}$文件定义的非线性向量函数$\text{c}(x),\text{ceq}(x)$；options定义了优化参数，可以使用$\text{MATLAB}$缺省参数设置。

例2 求下列非线性规划
$$\min f(x)=x_1^2+x_2^2+x_3^2+8$$

$$\text{s.t.}\{\begin{array}{l}\begin{array}{r}{x}_1^2-x_2+x_3^2\ge 0\\ x_1+x_2^2+x_3\le 20\\ -x_1-x_2+2=0\\ x_2+2x_3^2=3\\ x_1,x_2,x_3\ge 0\end{array}\end{array}$$

​ 解：

1. 编写$\text{M}$函数$\text{fun1.m}$定义目标函数：

function f=fun1(x)%定义目标函数
f=sum(x.^2)+8;
end

2. 编写$\text{M}$函数$\text{fun2.m}$定义非线性约束条件：

function [g,h]=fun2(x);%定义非线性约束条件
g=[-x(1)^2+x(2)-x(3)^2%回车换行，第一个不等式约束
    x(1)+x(2)^2+x(3)^3-20];%第二个不等式约束
h=[-x(1)-x(2)^2+2%回车换行，第一个等式约束
    x(2)+2*x(3)^2-3];%非线性等式约束，第二个等式约束
end

3. 编写主程序文件如下：

[x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[],'fun2');

​ 求得当$x_1=0.5522,x_2=1.2033,x_3=0.9478$时，最小值$y=10.6511$。这里的初始值是用了一组随机数，但是此题初始值不论是什么都不影响最终的结果。

二次规划

​ 若某非线性规划的目标函数为自变量的二次函数，约束条件又全是线性的，就称这种规划为二次规划。

$\text{MATLAB}$中的数学模型

$$\min \frac{1}{2}{\mathbf{x}}^{\text{T}}\mathbf{H}\mathbf{x}+{\mathbf{f}}^{\text{T}}\mathbf{x}$$

$$\text{s.t.}\{\begin{array}{l}\begin{array}{r}\mathbf{Ax}\le \mathbf{b}\\ \mathbf{Aeq}\cdot \mathbf{x}=\mathbf{beq}\\ \mathbf{lb}\le \mathbf{x}\le \mathbf{ub}\end{array}\end{array}$$

​ 这里$\mathbf{H}$是实对称矩阵，$\mathbf{f},\mathbf{b},\mathbf{beq},\mathbf{lb},\mathbf{ub}$是列向量，$\mathbf{A},\mathbf{Aeq}$是相应维数的矩阵。

求解命令

[x,fval]=quadporg(H,f,A,b,Aeq,beq,lb,ub,x0,options);

​ 返回值$\mathbf{x}$是决策向量$\mathbf{x}$的值，返回值$\text{fval}$是目标函数在$\mathbf{x}$处的值。

例3 求解二次规划
$$\min f(x)=2x_1^2-4x_1{x}_2+4x_2^2-6x_1-3x_2$$

$$\text{s.t.}\{\begin{array}{l}\begin{array}{r}{x}_1+x_2\le 3\\ 4x_1+x_2\le 9\\ x_1,x_2\ge 0\end{array}\end{array}$$

​ 解：编写如下程序：

H=[4,-4;-4,8];
f=[-6;-3];
A=[1,1;4,1];
b=[3;9];
lb=zeros(2,1);
[x,fval]=quadprog(H,f,A,b,[],[],lb,[],rand(2,1));

​ 求得$x_1=1.9500,x_2=1.0500,\min f(x)=-11.0250$，本题的最优解也与初始值无关。

应用实例

例4 供应与选址

​ 某公司有6个建筑工地要开工，每个工地的位置（用平面坐标系$a,\text{b}$表示，距离单位：千米）及水泥日用量$\text{d}$由下表给出。目前有两个临时料场位于$\text{A}(5,1),\text{B}(2,7)$，日储量各有20吨。假设从料场到工地之间均有直线道路相连。

1. 试制定每天的供应计划，即从$\text{A},\text{B}$两料场分别向各工地运送多少吨水泥，使总的吨千米数最小。
2. 为了进一步减少吨千米数，打算舍弃两个临时料场，改建两个新的，日储量各为20吨，问应建在何处，节省的吨千米数有多大？

工地位置（a,b）及水泥日用量d

	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.25
d	3	5	4	7	6	11

符号说明

$X_{\text{ij}}$	从料场$\text{j}$向工地$\text{i}$的运送量
$(a_{\text{i}},{\text{b}}_{\text{i}})$	工地$\text{i}$的位置，$\text{i}=1,2,\cdots ,6$
${\text{d}}_{\text{i}}$	工地$\text{i}$水泥的日用量，$\text{i}=1,2,\cdots ,6$
$(x_{\text{j}},y_{\text{j}})$	料场$\text{j}$的位置，$\text{j}=1,2$
${\text{e}}_{\text{j}}$	料场$\text{j}$的日储量，$\text{j}=1,2$

建立模型

使用临时料场的情形

​ 使用两个临时料场$\text{A}(5,1),\text{B}(2,7)$。在各工地用量必须满足和各料场不超过日储量的条件下，使总的吨千米数最小，这是线性规划问题，线性规划模型为：
$$\min f=\sum _{\text{j}=1}^2\sum _{\text{i}=1}^{6}aa(\text{i},\text{j})X_{\text{ij}}$$

$$\text{s.t.}\{\begin{array}{ll}\sum _{\text{j}=1}^2{X}_{\text{ij}}={\text{d}}_{\text{i}}& \text{i}=1,2,\cdots ,6\\ \sum _{\text{i}=1}^6{X}_{\text{ij}}\le {\text{e}}_{\text{i}}& \text{j}=1,2\end{array}$$

​ 其中$aa(\text{i},\text{j})=\sqrt{(x_{\text{j}}-a_{\text{i}})+(y_{\text{j}}-{\text{b}}_{\text{i}})}),\text{i}=1,2,\cdots ,6,\text{j}=1,2$，为常数。

不使用临时料场的情形

​ 由于料场的位置也是我们要求的决策变量，所以料场到各工地的距离不是常数了。这是非线性规划问题，但不是二次规划。
$$\min f=\sum _{\text{j}=1}^2\sum _{\text{i}=1}^6{X}_{\text{ij}}\sqrt{(x_{\text{j}}-a_{\text{i}})+(y_{\text{j}}-{\text{b}}_{\text{i}})})$$

$$\text{s.t.}\{\begin{array}{ll}\sum _{\text{j}=1}^2{X}_{\text{ij}}={\text{d}}_{\text{i}}& \text{i}=1,2,\cdots ,6\\ \sum _{\text{i}=1}^6{X}_{\text{ij}}\le {\text{e}}_{\text{i}}& \text{j}=1,2\end{array}$$

模型求解

使用临时料场的情形

​ 设:
$$\begin{gathered} X_{11}=X_1,X_{21}=X_2,X_{31}=X_3,X_{41}=X_4 \\ {X}_{51}=X_5,X_{61}=X_6,X_{12}=X_7,X_{22}=X_8 \\ {X}_{32}=X_9,X_{42}=X_{10},X_{52}=X_{11},X_{62}=X_{12} \end{gathered}$$
​ $\text{MATLAB}$程序如下：

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];%临时料场的日存储量

for i=1:6%计算每个工地到两个临时料场的距离
    for j=1:2
        aa(i,j)=sqrt((x(j)-a(i))^2+(y(j)-b(i))^2);
    end
end

c=[aa(:,1);aa(:,2)];%把aa矩阵拼接成列向量
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';
lb=zeros(12,1);
ub=[];
[x,fval]=linprog(c,A,B,Aeq,beq,lb,ub);

​ 计算结果为：

x=[3,5,0,7,0,1,0,0,4,0,6,10];
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

​ 总的吨千米数为136.2275。

不使用临时料场的情形

​ 设：
$$\begin{gathered} X_{11}=X_1,X_{21}=X_2,X_{31}=X_3,X_{41}=X_4 \\ {X}_{51}=X_5,X_{61}=X_6,X_{12}=X_7,X_{22}=X_8 \\ {X}_{32}=X_9,X_{42}=X_{10},X_{52}=X_{11},X_{62}=X_{12} \\ {x}_1=X_{13},y_1=X_{14},x_2=X_{15},y_2=X_{16} \end{gathered}$$
​ $\text{MATLAB}$程序如下：

1. 先编写M文件$\text{liaoch.m}$定义目标函数：

function f=liaoch(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];
f1=0;%料场1到各工厂的吨千米数
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;%料场2到各工厂的吨千米数
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

2. 编写源程序$\text{gying2.m}$：

d=[3 5 4 7 6 11];%工厂的水泥日用量
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=[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=d';
lb=[zeros(12,1);-inf;-inf;-inf;-inf];
ub=[];
[x,fval]=fmincon('liaoch',x0,A,B,Aeq,beq,lb,ub)

3. 计算结果为：

x=[2.9410 4.8405 3.8779 6.9431 1.3033 0.0221 0.0590 0.1595 0.1221 0.0569 4.6967   10.9779 5.7297 4.9757 7.2500 7.7500];
fval=90.4920;
exitflag=2;

​ 即两个新料场的坐标分别为$(5.7297,4.9757),(7.2500,7.7500)$，由料场A、B向6个工地运料方案为：

	1	2	3	4	5	6
料场A	2.9410	4.8405	3.8779	6.9431	1.3033	0.0221
料场B	0.0590	0.1595	0.1221	0.0569	4.6967	10.9779

​ 总的吨千米数为90.4920。比用临时料场节省了45.7355吨千米。

4. 若修改主程序$\text{gying2.m}$，取初值为上面的计算结果：

x0=[2.9410 4.8405 3.8779 6.9431 1.3033 0.0221 0.0590 0.1595 0.1221 0.0569 4.6967   10.9779 5.7297 4.9757 7.2500 7.7500];

​ 得计算结果如下：

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;

​ 此次计算结果与上次的结果不同，说明结果与初始值有关，初始值不同，则结果不同。并且此次的最优解要劣于上次的最优解。按照视频讲的是把一次结果的最优解作为初始值再求解一次，得到的最优解会变得更优，然后反复将最优解带入初始值，最优解会逐渐逼近唯一的最优解。但是这次最优解解要略差一些，个人认为这是因为精度的问题，因为每一次求解得到的决策变量都会保留很多位小数，而我带入的初始值舍弃了一部分小数，不过这两次的初始值得到的解相差不是很多，在误差的允许范围之内，可以接受。

数学建模老哥