python解决非线性优化问题 python非线性约束优化

软件库：scipy.optimize, numpy, CVXPY,Gekko 软件：octave 5.1，matlab

本文将介绍三种计算非线性约束优化的方法：
（1）scipy.optimize.minimize
（2）cvxpy
（3）Octave 5.1 sqp函数
（4）matlab ga函数

---

2020-08-13 更新

使用cvxpy库的时候，对于有些优化问题需要注意转换一下形式，例如下面二次规划问题(https://www.inverseproblem.co.nz/OPTI/index.php/Probs/QP)：

需要先将优化问题转成二次规划的一般形式：

$\text{min} f(\bold{x})=\frac{1}{2} \bold{x}$

其中$\bold{x}$是列向量，可以算出$\bold{H}=\begin{bmatrix} 1 & -1 \\ -1 & 2\end{bmatrix}$，$\bold{f}=\begin{bmatrix} -2 \\ -6 \end{bmatrix}$，相应的cvxpy计算程序如下（参考官方例子https://www.cvxpy.org/examples/basic/quadratic_program.html）：

import numpy as np
import cvxpy as cvx
x = cvx.Variable(2)
H = np.array([[1,-1],
             [-1,2]])
f = np.array([-2,-6])
A = np.array([[1,1],
     [-1,2],
     [2,1]])
b = np.array([2,2,3])

# obj = cvx.Minimize(1/2*x[0]**2+x[1]**2-x[0]*x[1]-2*x[0]-6*x[1])
obj = cvx.Minimize((1/2)*cvx.quad_form(x,H)+f.T@x)
cons = [A@x<=b, x[0]>=0]

prob = cvx.Problem(obj, cons)
prob.solve()

print(prob.status, prob.value, x.value)

在上面例子中，如果采用obj = cvx.Minimize(1/2*x[0]**2+x[1]**2-x[0]*x[1]-2*x[0]-6*x[1])计算会提示求解错误，需要处理为obj = cvx.Minimize((1/2)*cvx.quad_form(x,H)+f.T@x)形式，具体为什么前面的形式无法求解，笔者也还需进一步学习。

---

2019-10-19补充Gekko库非线性约束优化
使用Gekko优化计算库求解博文中的一个非线性约束优化例子（无法用cvxpy求解，报错提示优化问题不遵循DCP规则），程序代码如下：

from gekko import GEKKO
import time

start = time.time_ns()
m = GEKKO() # Initialize gekko
# Use IPOPT solver (default)
m.options.SOLVER = 3
# Change to parallel linear solver
# m.solver_options = ['linear_solver ma97']
# Initialize variables
x1 = m.Var(value=1,lb=0,ub=5)
x2 = m.Var(value=2,lb=0,ub=5)

# Equations
m.Equation(x1**2+x2**2<=6)
m.Equation(x1*x2==2)

m.Obj(x1**3-x2**2+x1*x2+2*x1**2) # Objective

m.options.IMODE = 3 # Steady state optimization


m.solve(disp=False) # Solve
print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
end = time.time_ns()

print('Objective: ' + str(m.options.objfcnval))
print('Time used: {:.2f} s'.format((end-start)*1e-9))
'''
Results
x1: [0.87403204792]
x2: [2.2882456138]
Objective: -1.040502879
Time used: 24.37 s
'''

相比于scipy.optimize.minimize和octave，使用gekko求解该问题耗时很久，估计是解法器没有选择到最优，因为刚接触gekko，该问题后面再研究。

Gekko默认使用公共服务器资源进行计算，网络传输会影响程序执行速度，在初始化Gekko时可设置为本地计算求解，只要设置m = Gekko(remote=False)即可，执行结果为：

---

MATLAB ga函数

2019-05-10更新

matlab有自带的遗传算法函数ga，相关语法如下图：

octave有个ga工具箱，同样提供了ga函数，程序和前面的类似，更改一下函数定义的位置即可。

对于下面octave中求解Himmelblau函数：
$f(x,y)=(x^{2} +y-11)^{2}+(x+y^2-7)^2$的最优解，在matlab中程序如下：

%% MATLAB ga求解最优化问题

lb=[-5,-5];
ub=[5,5]; 

%% initial guess
x0 = [-3,-2];
%% MATLAB ga函数，详细信息请查看帮助文档
x  = ga (@phi,2,[],[],[],[],lb,ub)

%% object function，matlab中函数定义要放电一个脚本的最后面，这和octave脚本函数位置要求不同
function obj = phi (x)
  obj = (x(1)^2+x(2)-11)^2+(x(1)+x(2)^2-7)^2;
end

结果为:

x = 1×2    
    3.0000    2.0000

从结构可以看出像ga这种进化算法只得到了一个最优解。因为matlab中的ga函数还可以处理线性约束和非线性约束问题，那么我们继续求解下文中的两个例子。

Example1

注意A，b系数矩阵要转换为Ax<=b的格式

% MATLAB ga求解最优化问题，注意A，b系数矩阵要转换为Ax<=b的格式
A = [-1,2;1,2;1,-2];
b = [2;6;2];
lb=[0,0];
ub=[]; 
option = optimoptions('ga','PlotFcn',@gaplotbestf)
%% MATLAB ga函数

x  = ga (@phi,2,A,b,[],[],lb,ub,[],option)

%% object function

function obj = phi (x)
  obj = (x(1)-1)^2+(x(2)-2.5)^2;
end

结果为：

x = 1×2    
    1.4224    1.7117

可见其结果离真实解还有一定偏差，为此我们可以控制ga的option来提供精度，将option修改为如下：

option = optimoptions('ga','ConstraintTolerance',1e-6,'PlotFcn',@gaplotbestf)

运行结果为：

x = 1×2    
    1.3941    1.6970

Example 2

option = optimoptions('ga','ConstraintTolerance',1e-6,'PlotFcn',@gaplotbestf)
%% MATLAB ga函数
x  = ga (@phi,2,[],[],[],[],[],[],@cons,option)

%% object function
function obj = phi (x)
  obj = x(1)^3-x(2)^3+x(1)*x(2)+2*x(1)^2;
end

%% nonlinear constraints，非线性约束条件
function [c,ceq] = cons(x)
    c = x(1)^2+x(2)^2-6; % Compute nonlinear inequalities at x，不等式约束，c(x)<=0
    ceq = x(1)*x(2)-2;       % Compute nonlinear equalities at x，等式约束，ceq(x)=0
end

结果如下：

Optimization terminated: stall generations limit exceeded

but constraints are not satisfied.

x =

1.0255    1.9512

可见计算结果并不满足约束条件，也和真实解（x1 = 0.8740320488968545
x2 = 2.2882456112715115）相差比较大。

ga作为一种随机性的进化算法，ga本身的参数（如crossover、mutation等）都会对结果造成很大影响，这里就不继续探讨了。

scipy.optimize.minimize

相关函数：
scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)

func：求解的目标函数，Python callable function
x0：给定的初始值
args：目标函数中的参数
method：求解算法
bounds：边界约束
constraints：条件约束

其他详细内容请查看官方帮助文档。

这部分内容为scipy官方文档翻译过来，但官方文档对带约束的最优化问题介绍的比较少，其引用了一本书籍的内容，下面记录之。

（1）根据文档所知，参考书籍为Nocedal, J, and S J Wright. 2006. Numerical Optimization. Springer New York. 使用例子来自example16.3(p462)，文档有误，问题描述如下：

代码如下：

import numpy as np
from scipy.optimize import minimize

# 目标函数
def objective(x):
    return (x[0] - 1)**2 + (x[1] - 2.5)**2

# 约束条件
def constraint1(x):
    return x[0] - 2 * x[1] + 2  #不等约束

def constraint2(x):
    
    return -x[0] - 2 * x[1] + 6 #不等约束

def constraint3(x):
        
    return -x[0] + 2 * x[1] + 2 #不等约束

# 初始猜想
n = 2
x0 = np.zeros(n)
x0[0] = 2
x0[1] = 0


# show initial objective
print('Initial SSE Objective: ' + str(objective(x0)))

# 边界约束
b = (0.0,None)
bnds = (b, b) # 注意是两个变量都要有边界约束

con1 = {'type': 'ineq', 'fun': constraint1} 
con2 = {'type': 'ineq', 'fun': constraint2}
con3 = {'type': 'ineq', 'fun': constraint3} 
cons = ([con1,con2,con3]) # 3个约束条件

# 优化计算
solution = minimize(objective,x0,method='SLSQP',\
                    bounds=bnds,constraints=cons)
x = solution.x

# show final objective
print('Final SSE Objective: ' + str(objective(x)))

# print solution
print('Solution')
print('x1 = ' + str(x[0]))
print('x2 = ' + str(x[1]))

结果如下：

Initial SSE Objective: 7.25
Final SSE Objective: 0.8000000011920985
Solution
x1 = 1.4000000033834283
x2 = 1.7000000009466527

就是当x1=1.4,x2=1.7时在约束条件下目标函数取得最小值。

再举一个例子加深印象，参照https://jingyan.baidu.com/album/6c67b1d69508b52787bb1edf.html?picindex=2，求解下面优化问题：

使用python代码如下：

import numpy as np
from scipy.optimize import minimize

def objective(x):
    return x[0]**3-x[1]**3+x[0]*x[1]+2*x[0]**2

def constraint1(x):
    return -x[0]**2 - x[1]**2 + 6.0

def constraint2(x):
    
    return 2.0-x[0]*x[1]

# initial guesses
n = 2
x0 = np.zeros(n)
x0[0] = 1
x0[1] = 2


# show initial objective
print('Initial SSE Objective: ' + str(objective(x0)))

# optimize
b = (0.0,3.0)
bnds = (b, b)
con1 = {'type': 'ineq', 'fun': constraint1} 
con2 = {'type': 'eq', 'fun': constraint2}
cons = ([con1,con2])
solution = minimize(objective,x0,method='SLSQP',\
                    bounds=bnds,constraints=cons)
x = solution.x

# show final objective
print('Final SSE Objective: ' + str(objective(x)))

# print solution
print('Solution')
print('x1 = ' + str(x[0]))
print('x2 = ' + str(x[1]))

结果如下：

Initial SSE Objective: -3.0
Final SSE Objective: -7.785844454002188
Solution
x1 = 0.8740320488968545
x2 = 2.2882456112715115

和原作者在matlab求解结果一致

特别注意： 使用约束的时候，要转化为形如ax+by>=c的形式，像前面第二个例子中的不等约束为$x^2+y^2<=6$，则要转化为$-x^2-y^2+6>=0$的形式，然后再编写约束函数的代码：

def constraint1(x):
    return -x[0]**2 - x[1]**2 + 6.0

接着定义约束类型：

con1 = {'type': 'ineq', 'fun': constraint1}

cvxpy

使用cvxpy求解前面第一个问题，起程序如下：

import cvxpy as cp
import numpy as np

'''
minimize f(x)=(x-1)^2+(y-2.5)^2
s.t.
x-2y+2>=0
-x-2y+6>=0
-x+2y+2>=0
x>=0
y>=0
'''
A = np.array([[1, -2],
            [-1, -2],
            [-1, 2]])

b = np.array([-2, -6, -2])

x = cp.Variable(2)
prob=cp.Problem(cp.Minimize( (x[0]-1)**2 + (x[1]-2.5)**2),
            [A@x>=b,x[0]>=0,x[1]>=0])
prob.solve()

# show final objective
print('Final SSE Objective: ' + str(prob.value))

# print solution
print('Solution of x and y: ' + str(x.value))

运行结果为：

Final SSE Objective: 0.8000000000000005
Solution of x and y: [1.4 1.7]

对于第二个例子笔者用cvxpy暂时解不出来，程序报错：

cvxpy.error.DCPError: Problem does not follow DCP rules. Specifically:
The objective is not DCP. Its following subexpressions are not:

因为求解不成功就不放出程序了，以后研究透了再更新。

---

Octave5.1

相关函数：

[…] = sqp (x0, phi, g, h, lb, ub, maxiter, tol)

该函数的官方文档介绍如下：

Minimize an objective function using sequential quadratic programming (SQP). 
Solve the nonlinear program 
min phi (x)
 x 
subject to 
g(x)  = 0
h(x) >= 0
lb <= x <= ub 
using a sequential quadratic programming method.

参数说明（重点）
（1）x0:初始猜想值，向量形式，如果是n个变量优化问题，x0也为n维向量；
（2）phi:目标函数，可以使用函数句柄
（3）g:等式约束条件，如果没有约束用空矩阵[]代替
（4）h：不等式约束条件，如果没有约束用空矩阵[]代替
（5）lb：变量左边界约束，同样注意和变量个数有关
（6）ub：变量右边界约束，和变量个数有关

先举个简单的例子，目标函数为Himmelblau函数，关于了解更多优化计算中测试算法使用的函数，可以参考网站http://infinity77.net/global_optimization/index.html下面的Test Function部分内容，也可在维基百科中搜索，这里就不过多介绍了。

Himmelblau函数：

$f(x,y)=(x^{2} +y-11)^{2}+(x+y^2-7)^2$

使用octave求解代码如下：

%%OCTAVE脚本里面有函数时，需要先有其他语句开始，不然就会当做是函数文件
1;
%% object function
function obj = phi (x)
  obj = (x(1)^2+x(2)-11)^2+(x(1)+x(2)^2-7)^2;
endfunction
lb=[-5,-5];
ub=[-5,5]; 

%% initial guess
x0 = [-3,-2];
%% 没有约束条件，用[]代替
[x, obj, info, iter, nf, lambda] = sqp (x0, @phi, [],[],lb,ub)

结果如下：

>> nonlinear_prob

x =

  -3.7793
  -3.2832

obj =    1.5554e-14
info =  104
iter =  8
nf =  17
lambda =

   0
   0
   0
   0

可见在初始值（-3，-2）情况下找到了最优解（-3.7793，-3.2832）

下面求解前面带约束的例子：

1;
%% equality constraints
function r = g (x)
  r = [ x(1)*x(2)-2];
endfunction

function r = h (x)
  r = [ 6-x(1)^2-x(2)^2];
endfunction

%% object function
function obj = phi (x)
  obj = x(1)^3-x(2)^3+x(1)*x(2)+2*x(1)^2;
endfunction
lb=[0,0];
ub=[3,3]; 

%% initial guess
x0 = [1,1];

[x, obj, info, iter, nf, lambda] = sqp (x0, @phi, @g, @h, lb, ub)

结果如下：

x =

   0.87403
   2.28825

obj = -7.7858
info =  104
iter =  6
nf =  8
lambda =

   7.03150
   4.58428
   0.00000
   0.00000
   0.00000
   0.00000

计算得出的x解和scipy.optimize求出的一致，这点以后有时间再验证一下。

————
2019-04-17补充
优路径最短问题可以归结为最优化问题，即求解一个点到其他点距离的最小值。

minimize

$minimize f(x,y)=\sum((x-datax_i)^2+(y-datay_i)^2)$
这里为了方便计算，使用了距离的平方，没有开根号计算。

根据前面的目标函数，我们写成符合scipy.optimzie.minimize的函数输入格式，最后完整的程序如下：

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

datax = 10*np.random.randn(15)
datay = 10*np.random.rand(15)

def f(x):
    return sum((x[0]-datax[i])**2+(x[1]-datay[i])**2 for i in range(len(datax)))

# 初值
x0 = [1,1]
# 计算结果
opt = minimize(f,x0=x0)
print(opt)

plt.figure()
plt.scatter(datax,datay,label='Original points')
plt.scatter(opt.x[0], opt.x[1], s=240, marker='h',label='Center point')
plt.legend()
# plt.show()
for i in range(len(datax)):
    plt.plot([opt.x[0],datax[i]],[opt.x[1],datay[i]])
plt.show()

结果如下：

fun: 762.1598697487602
 hess_inv: array([[ 3.33359580e-02, -8.50862577e-06],
       [-8.50862577e-06,  3.33328322e-02]])
      jac: array([ 0.00000000e+00, -7.62939453e-06])
  message: 'Optimization terminated successfully.'
     nfev: 28
      nit: 5
     njev: 7
   status: 0
  success: True
        x: array([0.89233268, 5.09306968])

CVXPY

import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt
datax = [2.5,6.7,3,9.3,9.2,9,4,7.2,0.3,2.8,10.3,2.4,5.7,8.1,3.5,6.89,1.34]
datay = [5.7,6.2,5.0,7.9,2.7,9.3,8,2,9,4,6.8,1.9,3.8,7.2,5.71,4.88,2.3]

def f(x):
    return sum((x[0]-datax[i])**2+(x[1]-datay[i])**2 for i in range(len(datax)))

x = cp.Variable(2)
prob=cp.Problem(cp.Minimize(sum((x[0]-datax[i])**2+(x[1]-datay[i])**2 for i in range(len(datax)))))
prob.solve()

# show final objective
print('Final SSE Objective: ' + str(prob.value))

# print solution
print('Solution of x and y: ' + str(x.value))

plt.figure()
plt.scatter(datax,datay,label='Original points')
plt.scatter(x.value[0], x.value[1], s=240, marker='h',label='Center point')
plt.legend()
# plt.show()
for i in range(len(datax)):
    plt.plot([x.value[0],datax[i]],[x.value[1],datay[i]])
plt.show()

对于求解线性规划问题，请参考文首给出的链接。

参考链接：
[1]https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#rdd2e1855725e-5
[2]https://apmonitor.com/pdc/index.php/Main/NonlinearProgramming
[3]https://jingyan.baidu.com/article/6c67b1d69508b52787bb1edf.html