Gekko/MATLAB/Scipy/CVXPY/Octave非线性约束优化

最新推荐文章于 2024-05-29 06:30:00 发布

落叶_小唱

最新推荐文章于 2024-05-29 06:30:00 发布

阅读量8.7k

点赞数 14

分类专栏： Python Math 文章标签：非线性优化 scipy.optimize minimize

本文链接：https://blog.csdn.net/ouening/article/details/83422337

版权

Python 同时被 2 个专栏收录

81 篇文章 14 订阅

订阅专栏

Math

37 篇文章 3 订阅

订阅专栏

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

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

博主其他相关博文：
CVXPY/Scipy/Octave线性规划问题求解

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{H}\bold{x}+\bold{f}'\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

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

落叶_小唱

关注

14
点赞
踩
47

收藏

觉得还不错? 一键收藏
4
评论
Gekko/MATLAB/Scipy/CVXPY/Octave非线性约束优化

软件库：scipy.optimize, numpy相关函数：scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)这部分内容为scipy官方...
复制链接

扫一扫