SciPy的optimize模块提供了许多数值优化算法，下面对其中的一些记录。
非线性规划（scipy.optimize.minimize）
一.背景：
现在项目上有一个用python 实现非线性规划的需求。非线性规划可以简单分两种，目标函数为凸函数 or 非凸函数。
凸函数的 非线性规划，比如fun=x2+y2+x*y，有很多常用的python库来完成，网上也有很多资料，比如CVXPY
非凸函数的 非线性规划（求极值），从处理方法来说，可以尝试以下几种：
1.纯数学方法，求导求极值；
2.使用神经网络，深度学习来处理，可参考反向传播算法中链式求导的过程；
3.寻找一些python库来做，本文介绍scipy.optimize.minimize的使用方法
二.库方法介绍
官方文档：https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
来看下改方法的入参
scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)
解释：
fun: 求最小值的目标函数
x0:变量的初始猜测值，如果有多个变量，需要给每个变量一个初始猜测值。minimize是局部最优的解法，所以
args:常数值,后面demo会讲解，fun中没有数字，都以变量的形式表示，对于常数项，需要在这里给值
method:求极值的方法，官方文档给了很多种。一般使用默认。每种方法我理解是计算误差，反向传播的方式不同而已，这块有很大理论研究空间
constraints:约束条件，针对fun中为参数的部分进行约束限制
三.demo
1.计算 1/x+x 的最小值
# coding=utf-8
from scipy.optimize import minimize
import numpy as np
#demo 1
#计算 1/x+x 的最小值
def fun(args):
a=args
v=lambda x:a/x[0] +x[0]
return v
if __name__ == "__main__":
args = (1) #a
x0 = np.asarray((2)) # 初始猜测值
res = minimize(fun(args), x0, method='SLSQP')
print(res.fun)
print(res.success)
print(res.x)
执行结果:函数的最小值为2点多，可以看出minimize求的局部最优
2.计算 (2+x1)/(1+x2) - 3x1+4x3 的最小值 x1,x2,x3的范围都在0.1到0.9 之间
# coding=utf-8
from scipy.optimize import minimize
import numpy as np
# demo 2
#计算 (2+x1)/(1+x2) - 3*x1+4*x3 的最小值 x1,x2,x3的范围都在0.1到0.9 之间
def fun(args):
a,b,c,d=args
v=lambda x: (a+x[0])/(b+x[1]) -c*x[0]+d*x[2]
return v
def con(args):
# 约束条件 分为eq 和ineq
#eq表示 函数结果等于0 ； ineq 表示 表达式大于等于0
x1min, x1max, x2min, x2max,x3min,x3max = args
cons = ({'type': 'ineq', 'fun': lambda x: x[0] - x1min},\
{'type': 'ineq', 'fun': lambda x: -x[0] + x1max},\
{'type': 'ineq', 'fun': lambda x: x[1] - x2min},\
{'type': 'ineq', 'fun': lambda x: -x[1] + x2max},\
{'type': 'ineq', 'fun': lambda x: x[2] - x3min},\
{'type': 'ineq', 'fun': lambda x: -x[2] + x3max})
return cons
if __name__ == "__main__":
#定义常量值
args = (2,1,3,4) #a,b,c,d
#设置参数范围/约束条件
args1 = (0.1,0.9,0.1, 0.9,0.1,0.9) #x1min, x1max, x2min, x2max
cons = con(args1)
#设置初始猜测值
x0 = np.asarray((0.5,0.5,0.5))
res = minimize(fun(args), x0, method='SLSQP',constraints=cons)
print(res.fun)
print(res.success)
print(res.x)
执行结果：
ADD:
全局最优的函数： scipy.optimize.basinhopping
有一个缺点是无法设置约束，求全局的最优解的函数
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.basinhopping.html
非线性方程组求解
SciPy中对非线性方程组求解是fslove()函数，它的调用形式一般为fslove(fun, x0)，fun是计算非线性方程组的误差函数，它需要一个参数x，fun依靠x来计算线性方程组的每个方程的值（或者叫误差），x0是x的一个初始值。
"""
计算非线性方程组：
5x1+3 = 0
4x0^2-2sin(x1x2)=0
x1x2-1.5=0
"""
## 误差函数
def fun(x):
x0,x1,x2 = x.tolist()
return[5*x1+3,4x0^2-2sin(x1x2),x1x2-1.5]
result = optimize.fsolve(fun,[1,1,1])
## result
[-0.70622057 -0.6 -2.5]
在计算非线性方程中的解时，比如像坐标上升算法，其中需要用到未知数的导数，同样，scipy的fslove()也提供了fprime参数传递未知数的雅各比矩阵从而加速计算，传递的雅各比矩阵每一行时某一方程对各个未知数的导数。对于上面的例子，我们可以写下如下的雅各比矩阵传入。
def j(x):
x0,x1,x2 = x.tolist()
return[[0,5,0],[8*x0,-2*x2*cos(x1*x2],[0,x2,x1]]
result = optimize.fsolve(fun,[1,1,1],fprime=j)
#result
[-0.70622057 -0.6 -2.5]
scipy的内部在实现fslove时应该时应该是利用了坐标上升算法或者梯度相关优化算法，但本人没有考证，有兴趣的可以看看源码。
最小二乘拟合
关于最小二乘算法的理论这里并不想谈，网上解释的文章也挺多，在 optimize模块中，可以使用leastsq()对数据进行最小二乘拟合计算。 leastsq()的用法很简单，只需要将计箅误差的函数和待确定参数的初始值传递给它即可。
x = np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])
y = np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])
def residual(p):
k,b = p
return y-(k*x+b)
r = optimize.leastsq(residual,[1,0])
k,b = r[0]
# print k
0.613495349193
# print b
1.79409254326
def func(x,p):
"""
计算的正弦波 ：A*sin(2*pi*k*x+theta)
"""
A,k,theta = p
return A*sin(2*np.pi*k*x+theta)
def redis(p,y,x):
return y-func(x,p)
x = np.linspace(0,2*np.pi,100)
A,k,theta = 10,0.34,np.pi/6
y0 = func(x,[A,k,theta])
# 加入噪声
np.random.seed(0)
y1 = y0+2*np.random.randn(len(x))
p0 = [7,0.40,0]
# p0是A,k,theta的初始值，y1，x要拟合的数据
plsq = optimize.leastsq(redis, p0,args=(y1,x))
print [A,k,theta] #真是的参数值
print plsq[0] #拟合后的参数值
对于像正弦波或者余弦波的曲线拟合，optimize提供curve_fit()函数，它的使用方式和leastq()稍有不同，它直接计算曲线的值，比如上面的拟合正弦波可以用cureve_fit()来写。
def func2(x,p):
"""
计算的正弦波 ：A*sin(2*pi*k*x+theta)
"""
A,k,theta = p
return A*sin(2*np.pi*k*x+theta)
ret,_=optimize.curve_fit(func2,x,y1,p0=p0)
该函数有一个缺点就是对于初始值敏感，如果初始频率和真实频率值差太多，会导致最后无法收敛到真是频率。
局部最小值
optimize模块还提供了常用的最小值算法如：Nelder-Mead、Powell、CG、BFGS、Newton-CG等，在这些最小值计算时，往往会传入一阶导数矩阵(雅各比矩阵)或者二阶导数矩阵(黑塞矩阵)从而加速收敛，这些最优化算法往往不能保证收敛到全局最小值，大部分会收敛到局部极小值。这些函数的调用方式为:
optimize.minimize(target_fun,init_val,method,jac,hess)
target_fun：函数的表达式计算；
init_val：初始值；
method：最小化的算法；
jac：雅各比矩阵
hess：黑塞矩阵。
全局最小值算法
全局最小值使用optimize.basinhopping()来实现，这个函数首先要定义一个误差计算方式，比如平方误差函数，niter时迭代的次数，最后还需要一个局部极小值优化方法，minimizer_kwargs传入。比如上面的正弦函数拟合：
def func1(x,p):
"""
计算的正弦波 ：A*sin(2*pi*k*x+theta)
"""
A,k,theta = p
return A*sin(2*np.pi*k*x+theta)
def func_error(p,y,x):
return np.sum((y-func1(x,p)**2)
result = optimize.basinhopping(func_error,[1,1,1],niter=10,
minimizer_kwargs={"method":"L-BFGS-B",
"args":(y1,x1)})
## [1,1,1]是传入的初始值，args是需要拟合的数据
题目：
1. 利用拉格朗日乘子法
#导入sympy包，用于求导，方程组求解等等
from sympy import *
#设置变量
x1 = symbols("x1")
x2 = symbols("x2")
alpha = symbols("alpha")
beta = symbols("beta")
#构造拉格朗日等式
L = 10 - x1*x1 - x2*x2 + alpha * (x1*x1 - x2) + beta * (x1 + x2)
#求导，构造KKT条件
difyL_x1 = diff(L, x1) #对变量x1求导
difyL_x2 = diff(L, x2) #对变量x2求导
difyL_beta = diff(L, beta) #对乘子beta求导
dualCpt = alpha * (x1 * x1 - x2) #对偶互补条件
#求解KKT等式
aa = solve([difyL_x1, difyL_x2, difyL_beta, dualCpt], [x1, x2, alpha, beta])
#打印结果，还需验证alpha>=0和不等式约束<=0
for i in aa:
if i[2] >= 0:
if (i[0]**2 - i[1]) <= 0:
print(i)
结果：
(-1, 1, 4, 6)
(0, 0, 0, 0)
2. scipy包里面的minimize函数求解
from scipy.optimize import minimize
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
#目标函数：
def func(args):
fun = lambda x: 10 - x[0]**2 - x[1]**2
return fun
#约束条件，包括等式约束和不等式约束
def con(args):
cons = ({'type': 'ineq', 'fun': lambda x: x[1]-x[0]**2},
{'type': 'eq', 'fun': lambda x: x[0]+x[1]})
return cons
#画三维模式图
def draw3D():
fig = plt.figure()
ax = Axes3D(fig)
x_arange = np.arange(-5.0, 5.0)
y_arange = np.arange(-5.0, 5.0)
X, Y = np.meshgrid(x_arange, y_arange)
Z1 = 10 - X**2 - Y**2
Z2 = Y - X**2
Z3 = X + Y
plt.xlabel('x')
plt.ylabel('y')
ax.plot_surface(X, Y, Z1, rstride=1, cstride=1, cmap='rainbow')
ax.plot_surface(X, Y, Z2, rstride=1, cstride=1, cmap='rainbow')
ax.plot_surface(X, Y, Z3, rstride=1, cstride=1, cmap='rainbow')
plt.show()
#画等高线图
def drawContour():
x_arange = np.linspace(-3.0, 4.0, 256)
y_arange = np.linspace(-3.0, 4.0, 256)
X, Y = np.meshgrid(x_arange, y_arange)
Z1 = 10 - X**2 - Y**2
Z2 = Y - X**2
Z3 = X + Y
plt.xlabel('x')
plt.ylabel('y')
plt.contourf(X, Y, Z1, 8, alpha=0.75, cmap='rainbow')
plt.contourf(X, Y, Z2, 8, alpha=0.75, cmap='rainbow')
plt.contourf(X, Y, Z3, 8, alpha=0.75, cmap='rainbow')
C1 = plt.contour(X, Y, Z1, 8, colors='black')
C2 = plt.contour(X, Y, Z2, 8, colors='blue')
C3 = plt.contour(X, Y, Z3, 8, colors='red')
plt.clabel(C1, inline=1, fontsize=10)
plt.clabel(C2, inline=1, fontsize=10)
plt.clabel(C3, inline=1, fontsize=10)
plt.show()
if __name__ == "__main__":
args = ()
args1 = ()
cons = con(args1)
x0 = np.array((1.0, 2.0)) #设置初始值，初始值的设置很重要，很容易收敛到另外的极值点中，建议多试几个值
#求解#
res = minimize(func(args), x0, method='SLSQP', constraints=cons)
#####
print(res.fun)
print(res.success)
print(res.x)
# draw3D()
drawContour()
结果：
7.99999990708696
True
[-1.00000002 1.00000002]
扫一扫
1674
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
