Python之建模规划篇--非线性规划

Posted 2022-12-10 风信子的猫Redamancy

Python之建模规划篇--非线性规划

• 基本介绍
• 线性规划与非线性规划的区别
• 非线性规划的Matlab解法
• Python 解决非线性规划
• • 1、等式约束下的拉格朗日乘子法
  • 2、Python实现对带约束的非线性规划求解
  • • Python编程实现求解
    • python使用SciPy库实现求解问题
    • 结果对比
  • 样例1
  • 样例2

基本介绍

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

对于一个实际问题，在把它归结成非线性规划问题时，一般要注意如下几点：
（i）确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基础上，确认什么是问题的可供选择的方案，并用一组变量来表示它们。
（ii）提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。
（iii）给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或“坏”的价值标准，并用某种数量形式来描述它。
（iv）寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果，因此还需要寻找出问题的所有限制条件，这些条件通常用变量之间的一些不等式或等式来表示。

线性规划与非线性规划的区别

如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到）；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。

非线性规划的Matlab解法

它的返回值是向量x，其中 FUN 是用M 文件定义的函数 f (x)；X0是x的初始值；A,B,Aeq,Beq 定义了线性约束 A* X ≤ B, Aeq * X = Beq ，如果没有线性约束，则A=[],B=[],Aeq=[],Beq=[]；LB 和UB 是变量x 的下界和上界，如果上界和下界没有约束，则LB=[]，UB=[]，如果x 无下界，则LB 的各分量都为-inf，如果x 无上界，则UB的各分量都为 inf；NONLCON 是用M 文件定义的非线性向量函数C(x),Ceq(x)；OPTIONS定义了优化参数，可以使用Matlab 缺省的参数设置。

下面给一个例子

多余的我就不多说了，matlab的可以详细看书
现在我们来看看再仔细看看非线性规划问题的解法到底有什么

Python 解决非线性规划

• 非线性规划可以简单分两种，目标函数为凸函数or非凸函数
• 凸函数的非线性规划，比如fun= x2 +y2 +xy,有很多常用库完成，比如cvxpy
• 非凸函数的非线性规划(求极值)，可以尝 试以下方法:
• • 纯数学方法，求导求极值
  • 神经网络、深度学习(反向传播算法中链式求导过程)
  • scipy. opt imize. minimize

scipy . optimize .minimize (fun, x0, args= () , method=None,jaC=None, hess=None, hes sp=None, bounds=None, constaints= ()，tol =None, callback=None, options=None)
fun:求最小值的目标函数
args:常数值
method:求极值方法，一般默认。
constraints:约束条件
x0:变量的初始猜测值，注意minimize是局部最优

1、等式约束下的拉格朗日乘子法

<font size=3公式推导：这里只是简单的放了一张图片对等式约束下拉格朗日乘子法的求解步骤进行了讲解。

2、Python实现对带约束的非线性规划求解

<font size=3求解实际例题

Python编程实现求解

#导入sympy包，用于求导，方程组求解等等
from sympy import *
#设置变量
x1 = symbols("x1")
x2 = symbols("x2")
alpha = symbols("alpha")
#beta = symbols("beta")

#构造拉格朗日等式
L = 60 - 10*x1 - 4*x2 + x1*x1 + x2*x2 - x1*x2 - alpha * (x1 + x2 - 8)

#求导，构造KKT条件
difyL_x1 = diff(L, x1) #对变量x1求导
difyL_x2 = diff(L, x2) #对变量x2求导
difyL_alpha = diff(L, alpha) #对alpha求导
#求解KKT等式
aa = solve([difyL_x1, difyL_x2, difyL_alpha], [x1, x2, alpha])
print(aa)

python使用SciPy库实现求解问题

from scipy.optimize import minimize
import numpy as np
#目标函数：
def func(args):
    fun = lambda x: 60 - 10*x[0] - 4*x[1] + x[0]**2 + x[1]**2 - x[0]*x[1] 
    return fun

#约束条件，包括等式约束和不等式约束
def con(args):
    cons = ('type': 'eq', 'fun': lambda x: x[0]+x[1]-8)
    return cons

if __name__ == "__main__":
    args = ()
    args1 = ()
    cons = con(args1)
    x0 = np.array((2.0, 1.0)) #设置初始值，初始值的设置很重要，很容易收敛到另外的极值点中，建议多试几个值

    #求解#
    res = minimize(func(args), x0, method='SLSQP', constraints=cons)
    print(res.fun)
    print(res.success)
    print(res.x)

结果对比

1.普通法求结果

2.scipy库求结果

通过两个Python程序的求解结果对比分析，两个的求解结果相差不大，并且误差都在一定范围内，所以可以认为两个的求解结果是一致的。

样例1

计算1/x+x的最小值
如果我们足够熟悉这个函数的性质，我们可以很容易得出来，在x>0的时候，他的最小值为2，因为这是一个双勾函数，为奇函数，可以由基本不等式得到他的最小值，但是我们怎么用程序去实现呢

from scipy.optimize import minimize
import numpy as np
#计算 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)
#     print(res.fun)
#     print(res.success)
#     print(res.x)

可以得到如下结果，我们可以得到函数的最小是为2点多，可以看出minimize求的局部最优

样例2

计算（2+x₁）/ （1+x₂) - 3x₁ + 4x₃的最小值，其中x₁、x₂、x₃范围在0.1到0.9之间

from scipy.optimize import minimize
import numpy as np

# 计算（2+x1）/ （1+x2) - 3x1 + 4x3的最小值，其中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)
    print(res.fun)
    print(res.success)
    print(res.x)

可以看出对于这类简单函数，局部最优解与真实最优解相差不大，但是对于复杂的函数，x0的初始值设置，会很大程度影响最优解的结果。

ADD:
全局最优的函数： scipy.optimize.basinhopping
有一个缺点是无法设置约束，求全局的最优解的函数
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.basinhopping.html

每日一句
Don’t try so hard, the best things come when you least expect them to.
(不要着急，最好的总会在最不经意的时候出现)