数模竞赛-非线性规划

这篇文章的Markdown我实在是修不好了…凑合着看吧。

基础概念

记$\overrightarrow x=[x_1,x_2,\cdots,x_n]^{\mathrm T}$是$n$维欧几里德空间$\mathbb{R}^n$中的一个点。$f\left(\overrightarrow x\right),g_i\left(\overrightarrow x\right),i=1,2,\cdots,p,h_j\left(\overrightarrow x\right),j=1,2,\cdots,q$是定义在$\mathbb{R}^n$上的实值函数。

若$f\left(\overrightarrow x\right),g_i\left(\overrightarrow x\right),i=1,2,\cdots,p,h_j\left(\overrightarrow x\right),j=1,2,\cdots,q$至少有一个是$\overrightarrow x$的非线性函数，则称如下形式的数学模型为非线性规划模型的一般形式：

$$
\displaylines{
    \min f\left(\overrightarrow x\right),\\
    \mathrm{s.t.}\begin{cases}
     g_i\left(\overrightarrow x\right)\leqslant0,&i=1,2,\cdots,p,\\
     h_j\left(\overrightarrow x\right)=0,&j=1,2,\cdots,q
    \end{cases}
}
$$

称满足所有约束条件的点$\overrightarrow x$的集合$K={x\in\mathbb{R}^n\mid g_i\left(\overrightarrow x\right)\leqslant0,i=1,\cdots,p;h_j\left(\overrightarrow x\right)=0,j=1,\cdots,q}$为非线性规划问题的约束集或可行域。$\forall\overrightarrow x\in K$称$\overrightarrow x$为非线性规划问题的可行解或可行点。

若$\overrightarrow{x}^\in K$，且$\forall\overrightarrow x\in K,f\left(\overrightarrow{x}^\right)\leqslant f\left(\overrightarrow x\right)$，则称$\overrightarrow{x}^*$为全局最优解，$f\left(\overrightarrow{x}^\right)$为全局最优值。若$\forall\overrightarrow x\in K,\overrightarrow x\neq\overrightarrow{x}^,f\left(\overrightarrow{x}^\right)<f\left(\overrightarrow x\right)$，则称$\overrightarrow{x}^$为严格全局最优解，$f\left(\overrightarrow{x}^*\right)$为严格全局最优值。

若 $\overrightarrow{x}^\in K$且存在$\overrightarrow{x}^$的邻域$N_{\delta}\left(\overrightarrow{x}^\right)$，$\forall\overrightarrow x\in N_{\delta}\left(\overrightarrow{x}^\right)\cap K,f\left(\overrightarrow{x}^\right)\leqslant f\left(\overrightarrow x\right)$，则称$\overrightarrow{x}^$为局部最优解，$f\left(\overrightarrow{x}^\right)$为局部最优值。若$\forall\overrightarrow x\in N_{\delta}\left(\overrightarrow{x}^\right)\cap K,\overrightarrow x\neq\overrightarrow{x}^*,f\left(\overrightarrow{x}^\right)<f\left(\overrightarrow x\right)$，则称$\overrightarrow{x}^$为严格局部最优解，$f\left(\overrightarrow{x}^*\right)$为严格局部最优解。

无约束非线性规划问题可具体表示为$\min f\left(\overrightarrow x\right),\overrightarrow x\in\mathbb{R}^n$。设$f\left(\overrightarrow x\right)$具有连续的一阶偏导数，$\overrightarrow{x}^*$是无约束问题的局部极小点，则$\nabla f\left(\overrightarrow{x}^*\right)=\overrightarrow0$。设矩阵$f\left(\overrightarrow x\right)$具有对各个变量的二阶偏导数，称下矩阵为$f\left(\overrightarrow x\right)$的Hesse矩阵，记为$\nabla^2f\left(\overrightarrow x\right)$。

$$
\left[
\begin{matrix}
\displaylines{
\dfrac{\partial^2f}{\partial x_1^2}&\dfrac{\partial^2f}{\partial x_1\partial x_2}&\cdots&\dfrac{\partial^2f}{\partial x_1\partial x_n}\\
\dfrac{\partial^2f}{\partial x_2\partial x_1}&\dfrac{\partial^2f}{\partial x_2^2}&\cdots&\dfrac{\partial^2f}{\partial x_2\partial x_n}\\
\vdots&\vdots&&\vdots\\
    \dfrac{\partial^2f}{\partial x_n\partial x_1}&\dfrac{\partial^2f}{\partial x_n\partial x_2}&\cdots&\dfrac{\partial^2f}{\partial x_n^2}
}
\end{matrix}
\right]
$$

设$f\left(\overrightarrow x\right)$有连续二阶偏导数，满足$\nabla f\left(\overrightarrow{x}^\right)=\overrightarrow 0$，且$\nabla^2f\left(\overrightarrow{x}^\right)$为正定阵，则$\overrightarrow{x}^*$为无约束优化问题的局部最优解。

设$\Omega$是$n$维欧几里德空间的一点集，若任意两点$\overrightarrow{x_1}\in\Omega,\overrightarrow{x_2}\in\Omega$，其连线上所有点$\alpha\overrightarrow{x_1}+(1-\alpha)\overrightarrow{x_2}\in\Omega,0\leqslant\alpha\leqslant1$，则称$\Omega$为凸集。两个凸集的交集是凸集。

给定函数$f\left(\overrightarrow x\right),\overrightarrow x\in D\subset\mathbb{R}^n$，若$\forall\overrightarrow{x_1},\overrightarrow{x_2}\in D,\lambda\in[0,1]$，有$f\left(\lambda\overrightarrow{x_1}+(1-\lambda)\overrightarrow{x_2}\right)\leqslant\lambda f\left(\overrightarrow{x_1}\right)+(1-\lambda)f\left(\overrightarrow{x_2}\right)$，则称$f\left(\overrightarrow x\right)$为$D$上的凸函数。

定义在凸集上的有限个凸函数的非负线性组合仍为凸函数。设$f\left(\overrightarrow x\right)$为定义在凸集$\Omega$上的凸函数，则对任一$\alpha\in\mathbb R$，集合$S_{\alpha}=\left{\overrightarrow x\mid\overrightarrow x\in\Omega,f\left(\overrightarrow x\right)\leqslant\alpha\right}$为凸集。

设$\Omega$为$n$维欧几里德空间$\mathbb{R}^n$上的开凸集，$f\left(\overrightarrow x\right)$在$\Omega$上具有一阶连续偏导数，则$f\left(\overrightarrow x\right)$为$\Omega$上的凸函数的充要条件是：$\forall\overrightarrow{x_1}\in\Omega,\overrightarrow{x_2}\in\Omega,\overrightarrow{x_1}\neq\overrightarrow{x_2},f\left(\overrightarrow{x_2}\right)\geqslant f\left(\overrightarrow{x_1}\right)+\nabla f\left(\overrightarrow{x_1}\right)^{\mathrm T}\left(\overrightarrow{x_2}-\overrightarrow{x_1}\right)$。若该式为严格不等式，则为严格凸函数的充要条件。

设$\Omega$为$n$维欧几里德空间$\mathbb{R}^n$上的开凸集，$f\left(\overrightarrow x\right)$在$\Omega$上具有二阶连续偏导数，则$f\left(\overrightarrow x\right)$为$\Omega$上的凸函数的充要条件是：$f\left(\overrightarrow x\right)$的Hesse矩阵$\nabla^2f\left(\overrightarrow x\right)$在$\Omega$上处处半正定。若$\forall x\in\Omega$，$f\left(\overrightarrow x\right)$的Hesse矩阵都是正定的，则$f\left(\overrightarrow x\right)$是$\Omega$上的严格凸函数。

若$f\left(\overrightarrow x\right)$为$\Omega$上的凸函数，$g_i\left(\overrightarrow x\right)$为$\mathbb{R}^n$上的凸函数，$h_j\left(\overrightarrow x\right)$为$\mathbb{R}^n$上的线性函数，则称该非线性规划问题为凸规划。

cvxpy只能求解凸规划。

非线性规划

例题：

求解非线性规划：

$$
\displaylines{
\min f\left(\overrightarrow x\right)=x_1^2+x_2^2-4x_1+4,\\
\mathrm{s.t.}\begin{cases}
g_1\left(\overrightarrow x\right)=-x_1+x_2-2\leqslant0,\\
g_2\left(\overrightarrow x\right)=x_1^2-x_2+1\leqslant0,\\
x_1,x_2\geqslant0
\end{cases}
}
$$

代码：

1
2
3
4
5
6
7
import numpy,cvxpy
x=cvxpy.Variable(2,pos=True)
obj=cvxpy.Minimize(sum(x**2)-4*x[0]+4)
con=[-x[0]+x[1]-2<=0,x[0]**2-x[1]+1<=0]
prob=cvxpy.Problem(obj,con)
prob.solve(solver='CVXOPT')
print(round(prob.value,4),numpy.round(x.value,4)) #最优值 最优解

二次规划

例题：

求解二次规划模型：

$$
\displaylines{
    \max -x_1^2-0.3x_1x_2-2x_2^2+98x_1+277x_2,\\
    \mathrm{s.t.}\begin{cases}
     x_1+x_2\leqslant100,\\
     x_1-2x_2\leqslant0,\\
     x_1,x_2\geqslant0
    \end{cases}
}
$$

求解，改为矩阵描述：

$$
\displaylines{
    \max\left[x_1,x_2\right]\left[
     \begin{matrix}
     -1&-0.15\\
     -0.15&-2
     \end{matrix}
    \right]\left[
     \begin{matrix}
     x_1\\x_2
     \end{matrix}
    \right]+\left[98,277\right]\left[
     \begin{matrix}
     x_1\\x_2
     \end{matrix}
    \right],\\
    \mathrm{s.t.}\begin{cases}
     \left[
     \begin{matrix}
     1&1\\
     1&-2
     \end{matrix}
     \right]\left[
     \begin{matrix}
     x_1\\x_2
     \end{matrix}
     \right]\leqslant\left[
     \begin{matrix}
     100\\0
     \end{matrix}
     \right],\\
     x_1,x_2\geqslant0
    \end{cases}
}
$$

代码：

1
2
3
4
5
6
7
8
9
10
11
import cvxpy,numpy
c2=numpy.array([[-1,-0.15],[-0.15,-2]])
c1=numpy.array([98,277])
a=numpy.array([[1,1],[1,-2]])
b=numpy.array([100,0])
x=cvxpy.Variable(2,pos=True)
obj=cvxpy.Maximize(cvxpy.quad_form(x,c2)+c1@x) #二次型
con=[a@x<=b]
prob=cvxpy.Problem(obj,con)
prob.solve(solver='CVXOPT')
print(numpy.round(x.value,4),round(prob.value,4)) #最优解 最优值

一般非线性规划求解

不是凸规划，就不能使用cvxpy求解。scipy.optimize.minimize函数需要满足这种标准型：

$$
\displaylines{
    \min f\left(\overrightarrow x\right),\\
    \mathrm{s.t.}\begin{cases}
     g_i\left(\overrightarrow x\right)\geqslant0,&i=1,2,\cdots,p,\\
     h_j\left(\overrightarrow x\right)=0,&j=1,2,\cdots,q
    \end{cases}
}
$$

例题1：

求解二次规划：

$$
\displaylines{
    \min -x_1^2-0.3x_1x_2-2x_2^2+98x_1+277x_2\\
    \mathrm{s.t.}\begin{cases}
     x_1+x_2\leqslant100,\\
     x_1-2x_2\leqslant0,\\
     x_1,x_2\geqslant0
    \end{cases}
}
$$

代码：

1
2
3
4
5
6
7
8
9
10
import numpy,scipy
c2=numpy.array([[-1,-0.15],[-0.15,-2]])
c1=numpy.array([98,277])
a=numpy.array([[1,1],[1,-2]])
b=numpy.array([100,0])
obj=lambda x:x@c2@x+c1@x
con={'type':'ineq','fun':lambda x:b-a@x} #>=不等式约束条件
bd=[(0,numpy.inf)for i in range(a.shape[1])] #取值范围
res=scipy.optimize.minimize(obj,numpy.random.randn(2),constraints=con,bounds=bd)
print(res) #fun是最优值 x是最优解

例题2：

求解二次规划：

$$
\displaylines{
    \min f\left(\overrightarrow x\right)=x_1^2+x_2^2+x_3^2+8,\\
    \mathrm{s.t.}\begin{cases}
     x_1^2-x_2+x_3^2\geqslant0,\\
     x_1+x_2^2+x_3^3\leqslant20,\\
     -x_1-x_2^2+2=0,\\
     x_2+2x_3^2=3,\\
     x_1,x_2,x_3\geqslant0
    \end{cases}
}
$$

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import numpy,scipy
obj=lambda x:sum(x**2)+8
def constr1(x):
    x1,x2,x3=x
    return [x1**2-x2+x3**3,20-x1-x2**2-x3**2]
def constr2(x):
    x1,x2,x3=x
    return [-x1-x2**2+2,x2+2*x3**2-3]
con1={'type':'ineq','fun':constr1}
con2={'type':'eq','fun':constr2} #等式约束条件
con=[con1,con2]
bd=[(0,numpy.inf)for i in range(3)]
res=scipy.optimize.minimize(obj,numpy.random.randn(3),constraints=con,bounds=bd)
print(res)

多目标规划

一般形式：

$$
\displaylines{
\min\overrightarrow f\left(\overrightarrow x\right)=\left[f_1\left(\overrightarrow x\right),f_2\left(\overrightarrow x\right),\cdots,f_m\left(\overrightarrow x\right)\right]^{\mathrm T},\\
\mathrm{s.t.}\begin{cases}
g_i\left(\overrightarrow x\right)\leqslant0,&i=1,2,\cdots,p,\\
h_j\left(\overrightarrow x\right)=0,& j=1,2,\cdots,q
\end{cases}
}
$$

其中$\overrightarrow x$为决策向量，$f_1\left(\overrightarrow x\right),f_2\left(\overrightarrow x\right),\cdots,f_m\left(\overrightarrow x\right)$为目标函数，记$\Omega=\left{\overrightarrow x\mid g_i\left(\overrightarrow x\right)\leqslant0,i=1,2,\cdots,p,h_j\left(\overrightarrow x\right)=0,j=1,2,\cdots,q\right}$为多目标规划的可行域或决策空间。$\overrightarrow f(\Omega)=\left{f\left(\overrightarrow x\right)\mid\overrightarrow x\in\Omega\right}$为多目标规划问题的像集或目标空间。

设$\overline{\overrightarrow x}\in\Omega$，若$\forall i=1,2,\cdots,m,\overrightarrow x\in\Omega;f_i\left(\overline{\overrightarrow x}\right)\leqslant f_i\left(\overrightarrow x\right)$。则称$\overline{\overrightarrow x}$为问题的绝对最优解，绝对最优解集为$\Omega_{ab}^*$。

设$\overline{\overrightarrow x}\in\Omega$，若不存在$\overrightarrow x\in\Omega$，使得$f_i\left(\overrightarrow x\right)\leqslant f_i\left(\overline{\overrightarrow x}\right),i=1,2,\cdots,m$，且至少有一个$f_j\left(\overrightarrow x\right)<f_j\left(\overline{\overrightarrow x}\right)$，则称$\overline{\overrightarrow x}$为问题的Pareto有效解，$f\left(\overline{\overrightarrow x}\right)$为有效点。

当$\overline{\overrightarrow x}\in\Omega$满足$f_i\left(\overline{\overrightarrow x}\right)\leqslant\alpha_i,i=1,2,\cdots,m$时，就认为$\overline{\overrightarrow x}$是一个满意解。

例题1：

$$
\displaylines{
    \min\left{-f_1\left(\overrightarrow x\right),f_2\left(\overrightarrow x\right)\right}\\
    \mathrm{s.t.}\begin{cases}
     0.5x_1+0.25x_2\leqslant8,\\
     0.2x_1+0.2x_2\leqslant4,\\
     x_1+5x_2\leqslant72,\\
     x_1+x_2\geqslant10,\\
     x_1,x_2\geqslant0
    \end{cases}
}
$$

其中利润极大化$\max f_1\left(\overrightarrow x\right)=2x_1+3x_2$，污染极小化$\min f_2\left(\overrightarrow x\right)=x_1+2x_2$。

法1：线性加权法

两个目标函数权重都取$0.5$，转化为：

$$
\displaylines{
    \min0.5\left(-2x_1-3x_2\right)+0.5\left(x_1+2x_2\right),\\
    \mathrm{s.t.}\begin{cases}
     0.5x_1+0.25x_2\leqslant8,\\
     0.2x_1+0.2x_2\leqslant4,\\
     x_1+5x_2\leqslant72,\\
     x_1+x_2\geqslant10,\\
     x_1,x_2\geqslant0
    \end{cases}
}
$$

法2：理想点法

先分别求俩公式各自的最优解，构造目标与这俩最优解的差的平方和。

先求这个MP的$f_1^*=-53$：

$$
\displaylines{
\min -2x_1-3x_2,\\
\mathrm{s.t.}\begin{cases}
0.5x_1+0.25x_2\leqslant8,\\
0.2x_1+0.2x_2\leqslant4,\\
x_1+5x_2\leqslant72,\\
x_1+x_2\geqslant10,\\
x_1,x_2\geqslant0
\end{cases}
}
$$

再求这个MP的$f_2^*=10$：

$$
\displaylines{
\min x_1+2x_2,\\
\mathrm{s.t.}\begin{cases}
     0.5x_1+0.25x_2\leqslant8,\\
     0.2x_1+0.2x_2\leqslant4,\\
     x_1+5x_2\leqslant72,\\
     x_1+x_2\geqslant10,\\
     x_1,x_2\geqslant0
    \end{cases}
}
$$

然后求这个：

$$
\displaylines{
    \min f=\left(-2x_1-3x_2+53\right)^2+\left(x_1+2x_2-10\right)^2,\\
    \mathrm{s.t.}\begin{cases}
     0.5x_1+0.25x_2\leqslant8,\\
     0.2x_1+0.2x_2\leqslant4,\\
     x_1+5x_2\leqslant72,\\
     x_1+x_2\geqslant10,\\
     x_1,x_2\geqslant0
    \end{cases}
}
$$

法3：序贯解法

当优先第一个目标函数时，即代入第一个目标函数并塞到第五个约束条件上：

$$
\displaylines{
    \min x_1+2x_2,\\
    \mathrm{s.t.}\begin{cases}
     0.5x_1+0.25x_2\leqslant8,\\
     0.2x_1+0.2x_2\leqslant4,\\
     x_1+5x_2\leqslant72,\\
     x_1+x_2\geqslant10,\\
     -2x_1-3x_1=-53,\\
     x_1,x_2\geqslant0
    \end{cases}
}
$$

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
import numpy,cvxpy

#法一
c1=numpy.array([-2,-3])
c2=numpy.array([1,2])
a=numpy.array([[0.5,0.25],[0.2,0.2],[1,5],[-1,-1]])
b=numpy.array([8,4,72,-10])
x=cvxpy.Variable(2,pos=True)
obj=cvxpy.Minimize(0.5*(c1+c2)@x)
con=[a@x<=b]
prob=cvxpy.Problem(obj,con)
prob.solve(solver='GLPK_MI')
print(x.value,prob.value) #最优解 最优值

#法2
obj1=cvxpy.Minimize(c1@x)
prob1=cvxpy.Problem(obj1,con)
prob1.solve(solver='GLPK_MI')
v1=prob1.value
obj2=cvxpy.Minimize(c2@x)
prob2=cvxpy.Problem(obj2,con)
prob2.solve(solver='GLPK_MI')
v2=prob2.value
print(v1,v2) #两个目标函数最优值
obj3=cvxpy.Minimize((c1@x-v1)**2+(c2@x-v2)**2)
prob3=cvxpy.Problem(obj3,con)
prob3.solve(solver='CVXOPT')
print(x.value) #最优解

#法三
con.append(c1@x==v1)
prob4=cvxpy.Problem(obj2,con)
prob4.solve(solver='GLPK_MI')
x3=x.value
print(x3) #最优解
print(-c1@x3) #利润
print(c2@x3) #污染物