线性规划的大M法和非线性规划的拉格朗日乘子法

一 、线性规划

什么是线性规划呢？
线性规划是运筹学的重要分支之一。(运筹学（operational research）是一门解决一定约束条件下最优解的学科，应用现有的科学技术知识与数学手段，来解决实际生活之中的各种问题，是一门应用学科。) 运筹学分支还有，规划论，排队论，图论，决策论等等。
在高中时我们已经接触过最简单的线性规划，如A,B两件商品的利润分别为2元与3元。同时满足商品A个数加商品B的个数不超过8个，A的个数不小于四个,B的个数不大于五个。我们可以设商品A的个数为x商品B的个数为 y 。所以上述问题可以表述为：目标函数： max Z=3 x1+2 x2约束条件：

解决这类问题的方法是根据约束条件画出可行域，一般可在可行域的顶点处取到最优解，这种方法叫做图解法。

链接：https://www.zhihu.com/question/320977814/answer/673901696

1.大M法

参考论文 用Excel演示大M单纯形法 楼建华 (石河子大学信息科学与技术学院*石河子832003)

大M单纯形法简称大M法。在大M法中，要求M足够大,通常,M作为符号参与运算。大M法单纯形中的数据均可表为aM + b的形式，如果用有序数对<a,b>等价表示aM +b,则大M法单纯形中的M被形式上消去,使得大M法可用Excel演示。

1.1 单纯形法

线性规划的标准形式为:
max CX

其中,B≥0,矩阵C、X、A、B的阶分别为1 x n、n x1，m x n,m x 1。
在有些教材中约定标准形为:求目标函数的最小值,A的秩r(A) = m,不要求B≥0。
单纯形的一s 般形式如表1所示。

其中,E≥0,(D E)是由(1)式中(A B)经初等行变换所得,(D E)中有m-r(A)行全为0(否则,(1)式没有可行解;邇常,删除单纯形中全为0的行),D中基变量所在列组成r(A)阶单位子矩阵,这个单位子矩阵称作基,同一基变量所在行列的交叉处为1通常,单纯形法由两个阶段组成:构造初始单纯形和换基。构造初始单纯形的一-般步骤为:对(1)式中(AB)做初等行变换，当变换所得(D E)中E≥0且D有r(A)阶单位子矩阵时，按表1相应填写初始单纯形中的XT、C、D、E、F、G、H.T,其中,T=C-FID。按上述步骤构造初始单纯形有一定难度,为此，对单纯形法做相应改进的算法应运而生，下述大M法就是其中之一-。

1.2 大M法

用单纯形法对下面的式子进行求解就是大M法

其中,Y=(y+2… y.)T是人工变量,N=(-M -M… -M),1是m阶单位矩阵, M充分大,其它部分与(1)成相同。
大M法依据(2)式和(1)式的下述关系,把求解(1)式的问题转换为求解(2)式:
①X0是(1)式最优解的充分必要条件是:(X0 0)是(2)式最优解。
②如果(2)式没有最优解,则(I)式也没有最优解。
③如果(2)式存在形如(X0 Y0)的最优解且Yo≠0,则(1)式没有最优解。
明显地,单位矩阵I可作为(2)式的初始基由它即可直接构造初始单纯形,简化了单纯形法的第一阶段,这是大M法的精华之所在。
对于(2)式,可采用表1形式的单纯形,但是,目标系数、基系数、检验数中含有M,而M很难预估，通常作为符号参与运算,无法用Excel 直接计算检验数和直接比较检验数的大小。
(2)式中,目标系数为两类,-类是C中的普通实数,另一类是N中的- M,均可表为aM + b的形式。检验数T=C-F"D,检验数均可表为aM+b的形式。以下用有序数对<a,b>等价表示aM +b,并采用表2所示的大M法单纯形:

其中,C、F、T分别是目标系数、基系数、检验数中M的系数,CO、F.T。分别是目标系数、基系数、检验数中所谓的常数项,T,=C-FTD,T.=Co-FD。形式上,表2中没有M,所有运算均为数值运算,因此,可用Excel演示大M法。
由于M充分大,当T，中最大数唯一时,无须计算T,当T,中最大数不唯- -时 ,只需计算T中最大数所在列的T。判断T,中最大数是否唯一的逻辑式比较繁琐，而且,运算量较大,以下改为:无论T,中最大数是否唯一,均计算T,中最大数所在列上的T0。
大M法的换基与普通单纯形法的换基相似,由四部分组成:检验最优解、选择人基变量、选择出基变量、换基。
最优解检验准则:当T1和T0中没有正数且基中人工变量均为0(所在行的约束值均为0)时,当前单纯形的可行基解是最优解(之一)。
人基准则:当T1或T0中有正数时,选择T。中最大数(不要求为正数)所在列入基,或者说，所在列的变量人基。
出基准则:当变量x人基时,如果x;的约束系数中没有正数,说明()式无上界,没有最优解。否则,x人基到{h, =e/d;ld;>0,i=1 ~m}中最小数所在行,该行上当前基变量出基。
换基:按照新基对应填写新的单纯形。当变量x人基到第i行时,C、F、F。的第i行相应改为x;和.x的目标系数,对(D E)做初等行变换,使D的第j列中第i行为1、其它行为0,计算新的检验数。

1.3用Excel演示大M法

例用大M法解下述线性规划:

第1步标准化。 引人松驰变量x、x以及人工变量y1 ,标准化为:

具体的运行步骤可以去参考这篇论文，这里直接展示结果。

用Excel自带的规划求解包求解

得到的结果和大M法手推一样，但是大M法太过于繁琐，一不小心就会犯错。

1.4用python实现大M法求解

直接使用python自带的库进行求解

代码如下

from scipy import optimize as op
import numpy as np
c=np.array([2,1,1])
A_ub=np.array([[0,-2,-1],[0,1,-1]])#不等式的系数
B_ub=np.array([2,1])#不等式的结果
A_eq=np.array([[1,-1,1]])#等式的系数
B_eq=np.array([2])#等式的结果
#确定范围
x1=(0,10)
x2=(0,10)
x3=(0,10)
res=op.linprog(c,A_ub,B_ub,A_eq,B_eq,bounds=(x1,x2,x3))
print(res)

结果如下

可以看到得到的结果和excle求解的几乎一致。
不使用库求解，这里是计算另一个方程式，我们只是为了比较一下python自带的库和自己编程用大M法求解的结果如何。

代码如下

import numpy as np

def pivot():
    l = list(d[0][:-2])
    jnum = l.index(max(l)) #转入编号
    m = []
    for i in range(bn):
        if d[i][jnum] == 0:
            m.append(0.)
        else:
            m.append(d[i][-1]/d[i][jnum])
    inum = m.index(min([x for x in m[1:] if x!=0]))  #转出下标
    s[inum-1] = jnum
    r = d[inum][jnum]
    d[inum] /= r
    for i in [x for x in range(bn) if x !=inum]:
        r = d[i][jnum]
        d[i] -= r * d[inum]
        
def solve():
    flag = True
    while flag:
        if max(list(d[0][:-1])) <= 0: #直至所有系数小于等于0
            flag = False
        else:
            pivot()
            
def printSol():
    for i in range(cn - 1):
        if i in s:
            print("x"+str(i)+"=%.2f" % d[s.index(i)+1][-1])
        else:
            print("x"+str(i)+"=0.00")
    print("objective is %.2f"%(-d[0][-1]))
d = np.loadtxt("大M法.txt", dtype=np.float)
(bn,cn) = d.shape
s = list(range(cn-bn,cn-1)) #基变量列表
solve()
printSol()

结果如下

使用库求解
代码如下

from scipy import optimize
import numpy as np
#确定c,A_ub,B_ub
c = np.array([1,14,6])
A_ub = np.array([[1,1,1],[1,0,0],[0,0,1],[0,3,1]])
B_ub = np.array([4,2,3,6])
res =optimize.linprog(-c,A_ub,B_ub)
print(res)

1.5 比较Excel手推，Excel自带包，python自带库，不使用python自带库四者大M法求解结果。

从上面的结果来看，四个方法得到的结果都相差无几，只是精度的问题，比较四个方式的操作难度，最难的肯定是用Excel自己推算大M法，要很细心才能得到正确的结果，用Excel自带的线性规划求解包做起来非常简单，一步到位。python两者肯定是使用自带的库的代码量要少一点，但是不用库更有利于自己理解算法。

二、非线性规划

非线性规划是具有非线性约束条件或目标函数的数学规划，是运筹学的一个重要分支。非线性规划是20世纪50年代才开始形成的一门新兴学科。70年代又得到进一步的发展。非线性规划在工程、管理、经济、科研、军事等方面都有广泛的应用，为最优设计提供了有力的工具。

1. 拉格朗日乘子法和KKT条件

拉格朗日乘数法是多元微分学中用来求函数z=f(x,y)在满足g(x,y)=0条件下的极值问题的方法：通过设F(x,y)=f(x,y)+λg(x,y),其中λ称为拉格朗日乘数，并求F(x,y)的极值点求得条件极值的方法
卡罗需－库恩－塔克条件
英文原名: Karush-Kuhn-Tucker Conditions常见别名:
Kuhn-Tucker，KKT条件，Karush-Kuhn-Tucker最优化条件，
Karush-Kuhn-Tucker条件，Kuhn-Tucker最优化条件，Kuhn-Tucker条件）是在满 足一些有规则的条件下，一个非线性规划（Nonlinear Programming）问题能有最优化解法的一个必要和充分条件。这是一个广义化拉格朗日乘数的成果。

KKT条件是指在满足一些有规则的条件下, 一个非线性规划(Nonlinear Programming)问题能有最优化解法的一个必要和充分条件. 这是一个广义化拉格朗日乘数的成果. 一般地, 一个最优化数学模型的列标准形式参考开头的式子, 所谓 Karush-Kuhn-Tucker 最优化条件，就是指上式的最优点x∗必须满足下面的条件:

1). 约束条件满足gi(x∗)≤0,i=1,2,…,p, 以及,hj(x∗)=0,j=1,2,…,q

2). ∇f(x∗)+∑i=1μi∇gi(x∗)+∑j=1λj∇hj(x∗)=0, 其中∇为梯度算子;

3). λj≠0且不等式约束条件满足μi≥0,μigi(x∗)=0,i=1,2,…,p。

KKT条件第一项是说最优点x∗必须满足所有等式及不等式限制条件, 也就是说最优点必须是一个可行解, 这一点自然是毋庸置疑的. 第二项表明在最优点x∗, ∇f必须是∇gi和∇hj的线性組合, μi和λj都叫作拉格朗日乘子. 所不同的是不等式限制条件有方向性, 所以每一个μi都必须大于或等于零, 而等式限制条件没有方向性，所以λj没有符号的限制, 其符号要视等式限制条件的写法而定.

为了更容易理解，我们先举一个例子来说明一下KKT条件的由来。

let L(x,μ)=f(x)+∑k=1μkgk(x)，其中μk≥0,gk(x)≤0

∵μk≥0 gk(x)≤0 => μg(x)≤0

∴maxμL(x,μ)=f(x) (2)

∴minxf(x)=minxmaxμL(x,μ) (3)

maxμminxL(x,μ)=maxμ[minxf(x)+minxμg(x)]=maxμminxf(x)+maxμminxμg(x)=minxf(x)+maxμminxμg(x)

又∵μk≥0, gk(x)≤0

　∴maxμminxμg(x)=0, 此时μ=0 or g(x)=0.

∴maxμminxL(x,μ)=minxf(x)+maxμminxμg(x)=minxf(x) (4)
　　此时μ=0 or g(x)=0.
　　联合(3),(4)我们得到minxmaxμL(x,μ)=maxμminxL(x,μ), 亦即

　　minxmaxμL(x,μ)=maxμminxL(x,μ)=minxf(x)

我们把maxμminxL(x,μ)称为原问题minxmaxμL(x,μ)的对偶问题，上式表明当满足一定条件时原问题、对偶的解、以及minxf(x)是相同的，且在最优解x∗处μ=0 or g(x∗)=0。把x∗代入(2)得maxμL(x∗,μ)=f(x∗)，由(4)得maxμminxL(x,μ)=f(x∗)，所以L(x∗,μ)=minxL(x,μ)，这说明x∗也是L(x,μ)的极值点，即
　　
参考链接： https://www.cnblogs.com/maybe2030/p/4946256.html
通常我们需要求解的最优化问题有如下几类：

(i) 无约束优化问题，可以写为:

                                  min f(x);  

(ii) 有等式约束的优化问题，可以写为:

                                   min f(x), 

                                        s.t. h_i(x) = 0; i =1, ..., n 

(iii) 有不等式约束的优化问题，可以写为：

                                  min f(x), 

                                        s.t. g_i(x) <= 0; i =1, ..., n

                                              h_j(x) = 0; j =1, ..., m

对于第(i)类的优化问题，常常使用的方法就是Fermat定理，即使用求取f(x)的导数，然后令其为零，可以求得候选最优值，再在这些候选值中验证；如果是凸函数，可以保证是最优解。

对于第(ii)类的优化问题，常常使用的方法就是拉格朗日乘子法（Lagrange Multiplier) ，即把等式约束h_i(x)用一个系数与f(x)写为一个式子，称为拉格朗日函数，而系数称为拉格朗日乘子。通过拉格朗日函数对各个变量求导，令其为零，可以求得候选值集合，然后验证求得最优值。

对于第(iii)类的优化问题，常常使用的方法就是KKT条件。同样地，我们把所有的等式、不等式约束与f(x)写为一个式子，也叫拉格朗日函数，系数也称拉格朗日乘子，通过一些条件，可以求出最优值的必要条件，这个条件称为KKT条件。

(a) 拉格朗日乘子法（Lagrange Multiplier)

对于等式约束，我们可以通过一个拉格朗日系数a 把等式约束和目标函数组合成为一个式子L(a, x) = f(x) + a*h(x), 这里把a和h(x)视为向量形式，a是横向量，h(x)为列向量，之所以这么写，完全是因为csdn很难写数学公式，只能将就了…。

然后求取最优值，可以通过对L(a,x)对各个参数求导取零，联立等式进行求取，这个在高等数学里面有讲，但是没有讲为什么这么做就可以，在后面，将简要介绍其思想。

(b) KKT条件

对于含有不等式约束的优化问题，如何求取最优值呢？常用的方法是KKT条件，同样地，把所有的不等式约束、等式约束和目标函数全部写为一个式子L(a, b, x)= f(x) + ag(x)+bh(x)，KKT条件是说最优值必须满足以下条件：

1. L(a, b, x)对x求导为零；

2. h(x) =0;

3. a*g(x) = 0;

求取这三个等式之后就能得到候选最优值。其中第三个式子非常有趣，因为g(x)<=0，如果要满足这个等式，必须a=0或者g(x)=0. 这是SVM的很多重要性质的来源，如支持向量的概念。

原文链接：https://blog.csdn.net/xianlingmao/article/details/7919597

2.拉格朗日乘子法计算方法

这里用一个例子来说明非等式约束如何计算的

第一步:加松弛变量x3、x4,使不等式约束变换为等式约束:

第二步:引入拉格朗日函数

第三步:引入新函数Z:

使用本程序时，先给定一个初始点X()，然后用计算机迭代计算，可求出最优解为:

再用一个等式约束例子来说明如何计算

计算过程

3.python求解拉格朗日乘子法

题目如下

代码如下

# 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)

运行结果如下：