文章目录
数学建模笔记（三）非线性规划和多目标优化
非线性规划概念和理论
非线性规划问题的数学模型
有约束非线性规划的求解
凸规划
一个简单的非线性规划模型
二次规划模型
多目标规划问题

数学建模笔记（三）非线性规划和多目标优化
非线性规划概念和理论
非线性规划问题的数学模型

有约束非线性规划的求解
常见的处理思路是：可能的话将非线性问题转化为线性问题，将约束问题转化为无约束问题。

1.对有等式约束的非线性规划问题使用Lagrange乘数法求解

2.罚函数法：

利用目标函数和约束函数构造带参数的增广目标函数，从而转换成一系列无约束非线性规划问题来进行求解。
凸规划
基本概念：

1.凸集与凸函数的定义

凸集：

凸函数：

例题1.

求解：

先对目标函数和约束函数中的非线性函数求二阶行列式，若均大于0，则是凸优化问题。

求解代码：
import numpy as np
import cvxpy as cp
x=cp.Variable(2,pos=True)
obj=cp.Minimize(sum(x**2)-4*x[0]+4)
con=[-x[0]+x[1]-2<=0,
x[0]**2-x[1]+1<=0]
prob = cp.Problem(obj, con)
prob.solve(solver='CVXOPT')
print("最优值为:",round(prob.value,4))
print("最优解为：\n", np.round(x.value,4))
#-------------结果--------------
最优值为: 3.7989
最优解为：
 [0.5536 1.3064]

一个简单的非线性规划模型
数学建模五步骤：

提出问题、选择建模方法、推导模型的数学表达式、求解模型、回答问题

例题2.

19英寸数量  $x_1$

21英寸数量  $x_2$

19英寸利润  $x_1* (339-195-0.01x_1-0.003x_2)$

21英寸利润  $x_2* (399-225-0.01x_2-0.004x_1)$

总利润  $x_1*(339-195-0.01x_1-0.003x_2)+x_2*(399-225-0.01x_2-0.004x_1)-400000$

问题求解：

import sympy as sp
import pylab as plt
plt.rc('text', usetex=True) #使用 LaTeX字体
plt.rc('font',size=14)
sp.var('x1, x2') #定义符号变量
y = (339-0.01*x1-0.003*x2)*x1+(399-0.004*x1-0.01*x2)*x2-(400000+195*x1+225*x2)
y = sp.simplify(y) #化简
dy1 = y.diff(x1) #求关于 x1的偏导
dy2 = y.diff(x2) #求关于 x2的偏导
s = sp.solve([dy1, dy2], [x1, x2])
x10 = round(float(s[x1])) #取整
x20 = round(float(s[x2]))
y0 = y.subs({x1: x10, x2: x20}) #符号函数代入数值
f = sp.lambdify('x1, x2', y, 'numpy') #符号函数转换为匿名函数
x = plt.linspace(0, 10000, 100)
X, Y = plt.meshgrid(x, x) #转换为网格数据
Z = f(X, Y)
ax=plt.subplot(121, projection='3d') #第一个子窗口三维画图
ax.plot_surface(X, Y, Z,cmap='viridis')
ax.set_xlabel('$x_1$'); ax.set_ylabel('$x_2$')
plt.subplot(122) #激活第二个子窗口
contr=plt.contour(X,Y,Z,10) #10条等高线
plt.clabel(contr) #等高线标注
plt.ylabel('$x_2$',rotation=0)
plt.xlabel('$x_1$')
sp.var('a', pos=True) #定义灵敏度 分析 的符号参数
y = (339-a*x1-0.003*x2)*x1+(399-0.004*x1-0.01*x2)*x2-(400000+195*x1+225*x2)
y = sp.simplify(y) #化简
dy1 = y.diff(x1) #求关于 x1的偏导
dy2 = y.diff(x2) #求关于 x2的偏导
s = sp.solve([dy1, dy2], [x1, x2])
sx1 = s[x1]; sx2 = s[x2] #提取解分量
s1 = sp.lambdify('a', sx1, 'numpy') #符号函数转换为匿名函数
s2 = sp.lambdify('a', sx2, 'numpy')
a0 = plt.linspace(0.002, 0.02, 50)
plt.figure()
plt.subplots_adjust(wspace = 0.65)
plt.subplot(121); plt.plot(a0, s1(a0))
plt.xlabel('$a$'); plt.ylabel('$x_1$')
plt.subplot(122); plt.plot(a0, s2(a0))
plt.xlabel('$a$'); plt.ylabel('$x_2$')
#plt.subplot_tool() #调整位置的工具
dx1 = sx1.diff(a); dx10 = dx1.subs(a, 0.01)
sx1a = dx10 * 0.01 / 4735
dx2 = sx2.diff(a); dx20 = dx2.subs(a, 0.01)
sx2a = dx20 * 0.01 / 7043
Y = y.subs({x1: s[x1], x2: s[x2]}) #求关于 a的目标函数
Y = sp.factor(Y); Y = sp.simplify(Y)
Ya = sp.lambdify('a', Y, 'numpy') #转换为匿名函数
a0 = plt.linspace(0.002, 0.02, 1000)
plt.figure(); plt.plot(a0, Ya(a0))
plt.xlabel('$a$'); plt.ylabel('$y$', rotation=0)
Sya = - 4735 ** 2 * 0.01 / 553641.025
y2 = y.subs({x1: 4735, x2: 7043, a: 0.011}) #计算近似最优利润
y3 = Y.subs(a, 0.011) #计算最优利润
delta = (y3 - y2) / y2 #计算利润的相对误差
plt.show()

二次规划模型
定义：

### 非线性规划的求解和应用

对于一般的非线性规划问题，由于不是凸优化，就不能使用cvxpy库求解。而需要使用scipy.optmize模块的minimize函数求解。

例题5.

由于不是凸函数，不能使用cvxpy库求解。

求解程序：
import numpy as np
from scipy.optimize import minimize
c2 = np.array([[-1, -0.15],[-0.15, -2]])
c1 = np.array([98, 277])
a = np.array([[1, 1], [1, -2]])
b = np.array([100, 0])
obj = lambda x: x @ c2 @ x + c1 @ x
con ={'type': 'ineq', 'fun': lambda x: b-a@x}
bd = [(0, np.inf) for i in range(a.shape[1])]
res = minimize(obj, np.ones(2), constraints=con, bounds=bd)
print(res) #输出解的信息
#-------------
 fun: 0.0
     jac: array([ 97.99999999, 276.99999997])
 message: 'Optimization terminated successfully'
    nfev: 6
     nit: 2
    njev: 2
  status: 0
 success: True
       x: array([0., 0.])

多目标规划问题
定义：

一般多目标规划问题没有最优解，只有满意解或有效解。

求有效解的预处理：

（1）无量纲化处理：每个目标函数的量纲通常是不一样的，在进行加权求解时由于量

纲的不可公度性，需要先进行无量纲化处理。

（2）数量级的归一化处理：当各个目标函数的数量级差异较大时，容易出现大数吃小

数现象，即数量级较大的目标在决策分析过程中容易占优，从而影响决策结果。

求解方法：

1.线性加权法

基本思想为根据目标的重要性确定一个权重，以目标函数的加权平均值为评价函数，是其达到最优。

2.$\epsilon$约束化

其中$\epsilon_i$是决策者对第i个目标容许接受阈值。

3.理想点法

基本思想：以每个单目标最优值为该目标的理想值，使每个目标函数值与理

想值的差的加权平方和最小。

4.优先级法

基本思想：根据目标重要性分成不同优先级，先求优先级高的目标函数的最优

值，在确保优先级高的目标获得不低于最优值的条件下，再求优先级低的目标函数。

    

                    

                    

                    
                    
                    
    

    与线性规划不同的是，非线性规划要求目标函数或约束条件中含有非线性函数。相应的求解这类问题就要用到非线性规划的方法。约束条件或者目标函数的放宽使得规划模型更具普适性，但也增加了问题求解的难度。对于简单的非线性规划问题，R语言中stat包即可求解。在这里我们给大家介绍R语言中求解非线性规划更为专业的Rdonlp2包。

    Rdonlp2在求解非线性规划问题上功能十分强大，用户可自行安装并查阅帮助文档：
library(Rdonlp2)
    Rdonlp2包的核心函数为 donlp2，可以快速求解非线性规划的最值。其用法如下：
donlp2(par,fn,
       par.upper=rep(+Inf,length(par)),
       par.lower=rep(+Inf,length(par)),
       A=NULL,
       lin.upper=rep(+Inf,length(par)),
       lin.lower=rep(-Inf,length(par)),
       nlin=list(),
       nlin.upper=rep(+Inf,length(nlin)),
       nlin.lower=rep(-Inf,length(nlin)),
       control=donlp2.control(),
       control.fun=function(lst){return(TRUE)},
       env=.GlobalEnv,
       name="Rdonlp2"
       )
    对该函数众多参数进行解释：
    par : 初始迭代向量
    fn : 连续型函数
    par.upper、par.lower : 自变量的上下限
    A：线性约束矩阵
    lin.upper、lin.lower : 线性约束条件的上下界
    nlin : 非线性约束条件函数列表
    nlin.upper、nlin.lower : 非线性约束条件的上下界
    其余参数可忽略。

    先看一例donlp2函数求解非线性规划的例子：   
    用Rdonlp2包编写非线性规划计算命令：
library(Rdonlp2)
p=c(0,0)              #迭代初始值
par.l=c(-100,100);par.u=c(100,100)#自变量定义域约束
fn=function(x){
   exp(x[1])*(4*x[1]^2+2*x[2]^2+4*x[1]*x[2]+2*x[2]+1)
}                  #目标函数
A=matrix(c(1,1,3,-1),2,byrow=TRUE)
lin.l=c(2,1);lin.u=c(+Inf,3)   #线性约束
nlcon1=function(x){
  -x[1]*x[2]
}                  #非线性约束
donlp2(p,fn,par.u,par.l,A,lin.l=lin.l,lin.u=lin.u,
       nlin=list(nlcon1))
       

$message
[1] "KT-conditions satisfied, no further correction computed"
$par
[1]  33.66667 100.00000
$gradf
[1]  1.625065e+19 2.243635e+17
    相应运行结果包括中间过程、算法参数、运行时间等数据，这里就不全部列出。
    再看一例利用Rdonlp2包求解二元rastrigin函数最小值的问题。该函数因为有众多极值点而仅有一个最值点对各类优化算法有一定的欺骗性，但用其来测试各类优化算法效果显著。该二元rastrigin函数为：
    先作出该函数的三维图像： 
a<-5
x<-seq(-a,a,0.01)
y<-seq(-a,a,0.01)
f<-function(x,y){
   x^2-10*sin(2*pi*x)+y^2-10*cos(2*pi*y)+20
}
z<-outer(x,y,f)
image(x,y,z,col=heat.colors(24))
library(rgl)
zorder<-rank(z)
persp3d(x,y,z,col=rainbow(as.integer(max(zorder)))[zorder])
fn<-function(x){
    f<-sum(x^2*cos(2*pi*x)+10)
}
par<-c(-500,500)
ret<-donlp2(par,fn)

ret$f
[1] 20
ret$par
[1] -2.654965e-06  2.654965e-06
    由输出结果可知该函数最小值为20 。所以Rdonlp2求解非线性规划也是非常方便的。当然了，还有很多更为一般的非线性规划问题是Rdonlp2包也不能解决的，这时候可能会向遗传算法、模拟退火这些启发式算法求助啦。这里且不作讨论。

    在许多实际问题中，衡量一个方案的好坏标准往往不止一个，例如制造一枚导弹，既要射程远又要最省燃料，还得精度高。这一类问题就不是单一的目标规划问题了，我们称之为多目标规划问题。
    求解多目标规划问题通常有主要目标法、分层序列法和线性加权求和法等方法。其中主要目标法通过确定一个主要目标，将多目标优化问题转化为线性或非线性规划问题；分层序列法是将目标中多个目标按照其重要程度排一个次序而逐步求解的过程；线性加权法则是对各目标赋予一个权数，加权求和得到一个新的目标函数而进行求解的过程。
    R语言中多目标规划问题的求解通常用goalprog包，其核心函数为llgp:
llgp(coefficient,targets,achievements,maxiter=1000,verbose=FALSE)
   其中coefficient为约束系数矩阵，targets为系数矩阵约束向量，achievements是包含objective、priority、p、n的目标函数，其中objective表示第几对偏差变量，priority表示该偏差变量的优先级，p和n分别为正负偏差变量的权系数。
    看一道多目标规划问题计算实例：（该题来自钱颂迪《运筹学》教材4.4节例5）
    llgp函数计算过程如下：
library(goalprog)
coefficients<-matrix(c(1,1,5,1,1,0,3,1),4)
targets<-c(10,4,56,12)
achievements<-data.frame(objective=1:4,priority=c(1,1,3,4),
                         p=c(2,3,0,1),n=c(0,0,1,0))
soln<-llgp(coefficients,targets,achievements)

soln$converged
[1] TRUE
soln$out
Decision variables
                X
X1   4.000000e+00
X2   6.000000e+00

Summary of objectives
        Objective           Over          Under         Target
G1   1.000000e+01   0.000000e+00   0.000000e+00   1.000000e+01
G2   4.000000e+00   0.000000e+00   0.000000e+00   4.000000e+00
G3   3.800000e+01   0.000000e+00   1.800000e+01   5.600000e+01
G4   1.000000e+01   0.000000e+00   2.000000e+00   1.200000e+01

Achievement function
                A
P1   0.000000e+00
P2   0.000000e+00
P3   1.800000e+01
P4   0.000000e+00
    由输出的结果我们可以看到x1、x2最优值分别为4和6。
    由该例我们可以看出，R语言求解多目标规划问题只需按照函数格式将规划的目标函数、约束条件套入格式即可，过程非常简单。

参考文献：
魏太云.R软件与最优化[C]//中国r语言会议.2008.





扫码关注数据科学家养成记






                
  

Python之建模规划篇--非线性规划
基本介绍
如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问
题。一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不象线性规划有
单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都
有自己特定的适用范围。
这是一个非线性规划问题的一般形式
对于一个实际问题，在把它归结成非线性规划问题时，一般要注意如下几点：
(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、等式约束下的拉格朗日乘子法
2、Python实现对带约束的非线性规划求解
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+x1)/ (1+x2) - 3x1 + 4x3的最小值，其中x1、x2、x3范围在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.
(不要着急，最好的总会在最不经意的时候出现)

整数规划

决策变量只能取整数的规划问题

动态规划

按时间或空间将多阶段决策问题划分为多个阶段，对每个阶段进行决策，前一次决策结果影响当前决策，当前决策影响下一次决策。



非线性规划

目标函数、约束条件至少一个是决策变量的非线性函数



可行解：满足约束条件的解

可行域：约束条件构成的区域


非线性规划问题包括无约束问题与有约束问题。

非线性规划问题引入了不等式约束，标志现代规划理论的开始。

无约束最优化

变量轮换法

最速下降法

牛顿法

共轭梯度法

变尺度法

有约束最优化

Zoutendijk可行方向法

近似规划法：在迭代点将非线性目标函数与约束条件一阶泰勒展开，转化为线性规划求解。求解的一系列解组成一个序列，序列收敛于原问题的解。

制约函数法：SUMT分为内点法（障碍函数）与外点法（惩罚函数）。内点法要求整个迭代都在可行域内进行，不能在边界上。外点法无要求，不在可行域都行。

二次规划QP：目标函数为决策变量的二次函数，约束函数是决策变量的线性函数

二次规划



某些问题就是求解一系列的二次规划问题，懂我再说什么问题，请留言。

正定二次规划

H为正定矩阵，是一个凸规划，局部最优解是全局最优解

一般正定二次规划问题


等式约束问题



用最优解方向将问题进一步转化



P为x出发到最优解x_的方向


所以目标函数写为


最终形式1


最终形式2


SQP逐步二次规划




将非线性规划转化为一系列的二次规划子问题求解。

B为H的近似矩阵，Bk+1通过BFGS公式由Bk修正得到

目标函数按照正定二次规划展开，约束条件取泰勒一阶近似