python绘制约束条件 python根据约束条件求解

转载

mob64ca14173efa 2023-10-30 20:03:01

文章标签 python绘制约束条件 optimization python 多项式优化问题 文章分类 Python 后端开发

Background

Analysis

1. python 符号运算

2. python嵌套优化问题

3. 将上面两部分结合起来

Reference

Background

信息论project，求LDPC的variable nodes和check nodes的degrees，也即λ和ρ，使得code rate最大，采用固定λ求ρ的方法（也可以固定ρ求λ，但要求max），抛开背景，就是个优化问题，形式如下：

$\quad \min \limits_{\rho _i} \sum_i \rho _i/ i \\ s.t. \left\{\begin{matrix}1 - \epsilon \lambda (1-\rho (x))\geq x,\forall x\in (1-\epsilon ,1) \\ \sum_i \rho_i=1 \end{matrix}\right.$

$where \left\{\begin{matrix} \sum_{i}\lambda _i = 1 \\ \rho(x)=\sum_i \rho_i x^{i-1} \\\lambda(x)=\sum_i \lambda_i x^{i-1} \end{matrix}\right.$

$\epsilon$

是给定常数，假设

$\lambda_i$

是自己设定好的（这里设置为

$\lambda_3=1, other\; \lambda_i=0$

, 即

$\lambda(x)=x^2$

），只需要求

$\rho_i$

Analysis

这个优化问题的第一条constraint非常特别，首先以一个多项式带入另一个多项式的x，自己算是没法算的，阶数非常高，因此要先用符号运算求出constraint的表达式，并且要想在一个区间内都满足，必须转换成在这个区间内的最小值大于等于0的一个求最小值的优化问题；

其次它在x的连续的区间内都要满足，引入了新的变量x，容易想到嵌套优化问题的思路，但直接scipy.optimize.minimize这个库会因为额外变量x而报错（具体原因我就不写了，试过就知道）， google了很久都没搜到，最后想出了一个解决方案，用优化问题定义另一个优化问题以作为前者的constraint。

1. python 符号运算

sympy库用来符号运算，定义好符号后，用这些符号作为多项式的系数，整个多项式就可以做运算了。这里顺便介绍一次定义多个符号的方法。

这里用scipy.poly1d 生成多项式，系数为p1,p2,p3,p4（scipy.poly1d就是numpy.poly1d），但用poly1d库生成的多项式虽然可以带入另一个多项式的x中，但不能print出表达式，看不到运算后的结果，因此这里选择用符号加法运算自己定义多项式。

from sympy import symbols
from scipy import poly1d

# 定义符号
x, p1, p2, p3 = symbols('x p1 p2 p3')

# 一次定义多个符号,返回tuple
vrs = symbols('p:5')  # (p0, p1, p2, p3, p4)
vrs = symbols('p1:5')  #    (p1, p2, p3, p4)

# 循环定义多组符号
for k in range(2,12):  
    multi_vrs = symbols('p1:'+ str(k))

# 1)用poly1d生成多项式，可拥有多项式任何属性，如获取系数、获取x为某个值时的结果等，
# 虽然可以带入另一个多项式的x，但不能print(G)
# ρ = ploy1d(list(vrs))
# print(ρ)  # p0*x**4 + p1*x**3 + p2*x**2 + p3*x + p4
# print(ρ.c)  # 获取系数
# print(ρ(1))  # x = 1时多项式的系数

# 2)自己定义多项式， sum运用了符号的加法运算来生成多项式
def polyn(x,vrs):
    n = len(vrs)
    result = 0 + sum([v * x**i for i, v in enumerate(vrs)])
    return result 
e = 0.4  # constant
l = [0, 0, 0, 1, 0, 0] # 假设λ多项式的系数为l3 = 1, 其他系数为0
λ = poly1d(l) # λ=x**2 (l3 = 1)

ρ = polyn(x, vrs)

G = 1 - e * λ(1 - ρ) - x
print(G)  # 1 - 0.4*(-p1 - p2*x - p3*x**2 - p4*x**3 + 1)**2

2. python嵌套优化问题

scipy.optimize库中，求解多元优化问题，可用fmin；如果带定义域（局部最小值），可用fminbound，只需设定x的取值区间；如果还要有constraints，就要用minimize。

这三个方法的优化目标函数都只能输入一个list作为唯一输入，这个list包含多个变量，即目标函数f(x), x[0], x[1], x[2]等都是变量，因此可以是多元的。

但是，如果一个优化问题的约束条件也是一个优化问题（如求某个函数的最小值，这个函数由原来优化的变量来确定），一定不要用fmin/fminbound这些封装好的库，因为没法重写源代码。

可以考虑用一个优化问题定义另一个优化问题，注意内层优化问题的返回值要取优化目标函数的值res.fun。同时还要注意求最大值的话要先给目标函数加上负号转为求最小值。

minimize的优化方法的选择有很多，分为Unconstrained minimization、Bound-Constrained minimization、Constrained Minimization等5种优化问题，不同类型的问题用的优化方法和需要的参数不同，可以根据官方文档的介绍选取。这里选取应用最广泛的'SLSQP'方法, 其他方法如TNC（Newton Conjugate-Gradient）也都可以试一下。

import numpy as np
from scipy.optimize import fminbound, minimize

# objective function: f
def f(l):
    n = len(l)
    l = [l[i] / (i + 1) for i in range(n)]  # pi / i, i = 1...n
    return sum(l)  # 注意加上负号转换为求最小值！！！
    
# inner_cons
def cons(l):
    # 用outer_minimize的变量l作为系数定义inner_minimize
    # 这里的多项式就是对应上面符号运算的结果，只把p换成l
    # G=1 - 0.4*(-p0*x**4 - p1*x**3 - p2*x**2 - p3*x - p4 + 1)**2
    funbnd = lambda x: 1 - 0.4*(-l[0] - l[1]*x - l[2]*x**2 - l[3]*x**3 + 1)**2 - x

    bnds2 = ((1 - e, 1),)
        
    # optimize f under cons
    res = minimize(funbnd, x0=(1 - e,), method='SLSQP', bounds=bnds2)
    # print("cons2:", res.fun)
    # print(res)
    return res.fun
   
# outer_cons
#* ineq means to be non-negative, eq means to be zero
Cons = (
    {'type': 'ineq', 'fun': cons}, # inner_cons
    
    {'type': 'eq', 'fun': lambda l:  sum(l) - 1}  # l1 + l2 + ... + ln = 1
)
    
# initial value of variables：array, tuple or list
initial_x = np.array([0 for _ in range(4)])
# initial_x = (0, 0, 0, 0)
# initial_x = [0, 0, 0, 0]

    
# bounds： array, tuple or List
bnds = ((0, 1), (0, 1),(0, 1), (0,1))

# outer_ominimize f under Cons
res = minimize(f, x0=initial_x, method='SLSQP', bounds=bnds, constraints=Cons)
print(res)

result: （2021.11.17:已修改为正确结果）

按照定义，objective function: f(l), l = [p1, p2, p3, p4].

nfev:30是迭代次数；fun:0.25是目标函数f(l)的最小值

x :array([0,0,0,1])是优化变量

$\rho_i$

的最终取值，分别对应p1=p2=p3=0,p4=1,对应的多项式是

$\rho(x)=0 x^0 + 0 x^1+ 0 x^2 + 1 x^3$

python绘制约束条件 python根据约束条件求解_python_10

LDPC(3,6) code, λ3=1时ρ6=1取得rate最大值0.5

python绘制约束条件 python根据约束条件求解_python_11

3. 将上面两部分结合起来

(1) 用sympy库中的lambdify将多项式定义为func(x)：

from sympy import symbols, lambdify
l = [0, 0, 1, 0, 0, 0]  # 多项式系数l : l1 x**0 + l2 x**1 +...+ln x**(n-1)
λ = polyn(x, l)  # 生成多项式 x**2
λ_func = lambdify(x, λ, 'numpy')  # 转为函数λ_func = x**2

(2)用symbol.subs()将多项式中的符号自动替换为某个特定的数值，这样就可以用循环全自动解多个优化问题了：

# inner_cons
    def cons(l):
        ## .subs([(p1,l[0], (p2,l[1])...(pn,l[n - 1])
        new_G = G.subs([(p[i], l[i]) for i in range(var_n)])  
        ## 用lambdify将多项式定义为函数，作为优化目标函数
        f = lambdify(x, new_G, 'numpy')

        bnds2 = ((1 - e, 1),)
        # inner optimize f under inner cons
        res = minimize(f, x0=(1 - e,), method='SLSQP', bounds=bnds2)

        return res.fun

2021.11.16 更新：

1）代码中优化目标函数的定义有小错误，但不影响变量p1,p2,p3,p4的取值；

2）代码中重新定义(p1,p2,p3...pn)为(x^0,x^1,x^2...x^n)的系数，便于理清思路；

3）注意求最大值要给目标函数加负号转为求最小值;

4) 代码中约束条件G的表达式有修改。

2021.11.17 更新：

1）将前两个部分结合起来。

2) 踩坑：我也是第一次用minimize库，这里给大家表演了一个踩坑：

minimize库只能求local minimal，它的初始值x0直接影响求得的结果，因为它是基于这个初始点沿着梯度方向去搜索的；求全局最优解的库叫scipy.basinhopping，但它不能解带约束条件的优化问题，而且如果变量有bounds，遇到复杂的函数也经常解不出来。

因此如果要解带约束的问题，可以随机生成多个初始点分布求最优解，再取最优值。

另外，minimize有的优化问题也是解不出来的，要看message是不是'Optimization terminated successfully'。

2021.11.22 更新：

固定λ优化ρ的完整代码：

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import fminbound, minimize
from scipy import poly1d
from sympy import symbols, lambdify

np.set_printoptions(precision=3, suppress = True,formatter={'float': '{: 0.3f}'.format}) # turn off np’s scientific notation
    
e = 0.4
x = symbols('x')


def polyn(x,vrs,length):
    # order : p1 x^1 + p2 x^2 + p3 x^3 + pn x^n
    result = 0 + sum([vrs[i] * x**i for i in range(length)])
    return result 
    

def find_min(G, p, l, var_n):
    # inner_cons: method1
    # def cons(l):
    #     new_G = G.subs([(p[i], l[i]) for i in range(var_n)])
    #     f = lambdify(x, new_G, 'numpy')
    #     bnds2 = ((1 - e, 1),)
    #     # inner optimize f under inner cons
    #     res1 = minimize(f, x0=((2 - e) / 2,), method='SLSQP', bounds=bnds2)
    #     res2 = minimize(f, x0=(1 - e,), method='SLSQP', bounds=bnds2)
    #     res3 = minimize(f, x0=(1,), method='SLSQP', bounds=bnds2)
    #     # print(res.fun, res.x)
    #     # print(res.message)
    #     return min(res1.fun, res2.fun, res3.fun)
    
    # inner_cons: method2
    def cons(l):
        new_G = G.subs([(p[i], l[i]) for i in range(var_n)])
        f = lambdify(x, new_G, 'numpy')
        sample_point = np.linspace(1-e,1,num=1000)
        # min = 0
        y = f(sample_point)
        return min(y)
   
    # outer_cons
    #* ineq means to be non-negative, eq means to be zero
    Cons = (
        {'type': 'ineq', 'fun': cons}, # G_min >= 0
        
        {'type': 'eq', 'fun': lambda l:  sum(l) - 1}  # l1 + l2 + ... + ln = 1
    )
    
    # objective function: f
    def f(l):
        return sum([l[i] / (i + 1) for i in range(len(l))])
        
    # initial guess
    initial_x = np.array([1 / var_n for _ in range(var_n)])
    
    # bounds
    bnds = [(0, 1) for _ in range(var_n)]

    # outer optimize f under outer Cons
    # f(l), l = [l[0], l[1], l[2], l[3],...]
    res = minimize(f, x0=initial_x, method='SLSQP', bounds=bnds, constraints=Cons) #  SLSQP, Newton-CG, 'BFGS'
    print(res)

    sum_ρ = res.fun
    sum_λ = f(l)
    rate = 1 - sum_ρ / sum_λ
    print(rate)

    return rate


if __name__ == '__main__':
    # λ(x)
    # l = [0.000,  0.417,  0.068,  0.117,  0.268,  0.108,  0.022]
    l = [0, 0, 1, 0, 0, 0]
    λ = polyn(x, l, len(l))
    λ_func = lambdify(x, λ, 'numpy')
    
    total_p_len = 20
    p = symbols(f"p1:{total_p_len + 1}")
 
    max_rate = 0
    j = 0
    for i in range(7, 13):  # len = 7~12
        ρ = polyn(x, p, i)
        # print(ρ)
        # constraint expression
        G = 1 - e * λ_func(1 - ρ) - x
        print(G)
        res = find_min(G, p, l, i)
        if res > max_rate:
            max_rate = res
            j = i
    print(j, max_rate)