给书配代码-电力经济调度（1）：基于拉格朗日及运筹规划方法的经济调度算法

小咪咪 2022-10-30 15:24 231阅读 0赞

电气博文传送门  
[学好电气全靠它，个人电气博文目录（持续更新中…）][Link 1]

## 知识点 ##

电力经济调度(Economic Dispatch, ED)的目标是追求某个研究时段内所有开机机组(Committed Unit)的运行成本最小，约束条件包括电力供需平衡及机组出力限制，优化变量是每台机组的出力。由于ED一般是在机组组合(Unit Commitment, UC)结束之后进行，因此如下的经济模型中不会包含反映机组开、停状态的0/1二元整型变量。  
经济模型中个会包含反映机组开、停状态的0/1 二元整型变量。  
![2021022009212854.png_pic_center][]

![20210220092723638.png_pic_center][] ![20210220092802271.png_pic_center][] ![20210220093528729.png_pic_center][] ![20210220093528729.png_pic_center][] ![20210220093614722.png_pic_center][]

## 算例 ##

![20210220093804128.png_pic_center][] ![20210220093831278.png_pic_center][]

## 程序python实现 ##

非线性规划（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中为参数的部分进行约束限制

**常规解法**  
语言：python

from scipy.optimize import minimize
    import numpy as np
    #目标函数即min(FG1+FG2+FG3)
    def fun(args):
        a1,a2,a3,b1,b2,b3,c1,c2,c3=args
        v=lambda x: (a1+a2*x[0]+a3*x[0]*x[0]+b1+b2*x[1]+b3*x[1]*x[1]+c1+c2*x[2]+c3*x[2]*x[2])
        return v
    
    def con(args):
        # 约束条件 分为eq 和ineq
        #eq表示 函数结果等于0 ； ineq 表示 表达式大于等于0
        d, x0min, x0max,x1min, x1max,x2min,x2max= args1
        # (700, 100, 200, 120, 250, 150,300) # d1, x0min, x0max,x1min, x1max,x2min,x2max
    
        #约束条件d1-x[0]-x[1]-x[2]转换为x[0]+x[1]+x[2]=700
        #100<=x[0]<=200,120<=x[1]<=250,150<=x[2]<=300
    
        cons = ({ 'type': 'eq', 'fun': lambda x: d-x[0]-x[1]-x[2]},\
                  { 'type': 'ineq', 'fun': lambda x:x[0]-x0min}, \
                { 'type': 'ineq', 'fun': lambda x: -x[0] + x0max}, \
                { 'type': 'ineq', 'fun': lambda x: x[1] - x1min}, \
                { 'type': 'ineq', 'fun': lambda x: -x[1] + x1max}, \
                { 'type': 'ineq', 'fun': lambda x: x[2] - x2min}, \
                { 'type': 'ineq', 'fun': lambda x: -x[2] + x2max})
    
        return cons
    
    
    if __name__ == "__main__":
        # 定义常量值
        args = (4,0.3,0.0007,3,0.32,0.0004,3.5,0.3,0.00045)  # a1,a2,a3,b1,b2,b3,c1,c2,c3
        # 设置参数范围/约束条件
        args1 = (700, 100, 200, 120, 250, 150,300)  # d1, x0min, x0max,x1min, x1max,x2min,x2max
        cons = con(args1)
        # 设置x初始猜测值
        x0 = np.array((150, 250, 20))
        res = minimize(fun(args), x0, method='SLSQP', constraints=cons)
        print('代价',res.fun)
        print(res.success)
        print('解',res.x)

![20210220100525226.png_pic_center][]

如果使用书上的解得到的代价为4+0.3x175+0.0007x175x1 75+3+0.32x250+0.0004x250x250+3.5+0.3x275+0.00045*275*275=305.96875

**书上解法复现**  
1、求解参数

import matplotlib.pyplot  as plt
    from pylab import *
    mpl.rcParams['font.sans-serif'] = ['SimHei']
    mpl.rcParams['axes.unicode_minus'] = False
    import numpy as np
    
    x0=np.arange(100,201,1)#100-200的整数
    y=4+0.3*x0+0.0007*x0*x0
    plt.plot(x0,y)
    plt.title('第一台机组运行曲线FG1')
    plt.show()
    #线性规划参数
    print('b',4+0.3*100+0.0007*100*100)
    print('PI11阶段参数',((4+0.3*125+0.0007*125*125)-(4+0.3*100+0.0007*100*100))/(125-100))
    print('PI12阶段参数',((4+0.3*150+0.0007*150*150)-(4+0.3*125+0.0007*125*125))/(150-125))
    print('PI13阶段参数',((4+0.3*175+0.0007*175*175)-(4+0.3*150+0.0007*150*150))/(175-150))
    print('PI14阶段参数',((4+0.3*200+0.0007*200*200)-(4+0.3*175+0.0007*175*175))/(200-175))

![20210220105800208.png_pic_center][] 求解得到的参数和书上给的解很接近。 ![2021022010584397.png_pic_center][]

同理得到FG2,

x0=np.arange(120,251,1)#120-250的整数
    y=3+0.32*x0+0.0004*x0*x0
    plt.plot(x0,y)
    plt.title('第二台机组运行曲线FG2')
    plt.show()
    
    #线性规划参数
    print('b',3+0.32*120+0.0004*120*120)
    print('PI21阶段参数',((3+0.32*220+0.0004*220*220)-(4+0.32*100+0.0004*100*100))/(220-100))
    print('PI22阶段参数',((3+0.32*250+0.0004*250*250)-(3+0.32*220+0.0004*220*220))/(250-220))

![20210220110600690.png_pic_center][] ![20210220110632176.png_pic_center][]

FG3的参数

x0=np.arange(150,301,1)#150-300的整数
    y=3.5+0.3*x0+0.00045*x0*x0
    plt.plot(x0,y)
    plt.title('第三台机组运行曲线FG3')
    plt.show()
    
    #线性规划参数
    print('b',3.5+0.3*150+0.00045*150*150)
    print('PI31阶段参数',((3.5+0.3*200+0.00045*200*200)-(3.5+0.3*150+0.00045*150*150))/(200-150))
    print('PI32阶段参数',((3.5+0.3*250+0.00045*250*250)-(3.5+0.3*200+0.00045*200*200))/(250-200))
    print('PI33阶段参数',((3.5+0.3*300+0.00045*300*300)-(3.5+0.3*250+0.00045*250*250))/(300-250))

![20210220113624421.png_pic_center][] ![20210220113706280.png_pic_center][]

我们求得的运行参数和书上给的，只是在位数上不同。为了统一，我后面采取书上的参数。  
2、运行求解

from scipy.optimize import minimize
    import numpy as np
    #目标函数即min(FG1+FG2+FG3)
    def fun(args):
        Z0,a0,a1,a2,a3,Z1,b0,b1,Z2,c0,c1,c2=args
        #Z0+x[0]*a0+x[1]*a1+x[2]*a2+x[3]*a3为FG1
        #Z1+x[4]*b0+x[5]*b1为FG2
        #Z2+x[6]*c0+x[7]*c1+x[8]*c2 为FG3
        v=lambda x: (Z0+x[0]*a0+x[1]*a1+x[2]*a2+x[3]*a3+Z1+x[4]*b0+x[5]*b1+Z2+x[6]*c0+x[7]*c1+x[8]*c2)
        return v
    
    def con(args):
        # 约束条件 分为eq 和ineq
        #eq表示 函数结果等于0 ； ineq 表示 表达式大于等于0
        d, x1min, x1max,x21min, x21max,x22min,x22max,x3min,x3max= args1
    
        #等式约束条件d-x[0]-x[1]-x[2]-x[3]-x[4]-x[5]-x[6]-x[7]-x[8]=700-100-120-150-x[0]-x[1]-x[2]-x[3]-x[4]-x[5]-x[6]-x[7]-x[8]
        #不等式约束条件
        #0<=x[0,1,2,3]<=25 x1min=0,x1max=25 PG1的粤省事
        #0<=x[4]<=100 x21min=0,x21max=100,0<=x[5]<=30 x22min=0,x22max=30 PG2的约束
        #0<=x[6,7,8]<=50 x3min=0,x3max=50 PG3的约束
        cons = ({ 'type': 'eq', 'fun': lambda x,: d-x[0]-x[1]-x[2]-x[3]-x[4]-x[5]-x[6]-x[7]-x[8]},\
                  { 'type': 'ineq', 'fun': lambda x:x[0]-x1min},
                { 'type': 'ineq', 'fun': lambda x: -x[0] + x1max},
                { 'type': 'ineq', 'fun': lambda x: x[1] - x1min},
                { 'type': 'ineq', 'fun': lambda x: -x[1] + x1max},
                { 'type': 'ineq', 'fun': lambda x: x[2] - x1min},
                { 'type': 'ineq', 'fun': lambda x: -x[2] + x1max},
                { 'type': 'ineq', 'fun': lambda x: x[3] - x1min},
                { 'type': 'ineq', 'fun': lambda x: -x[3] + x1max},
                { 'type': 'ineq', 'fun': lambda x: x[4] - x21min},
                { 'type': 'ineq', 'fun': lambda x: -x[4] + x21max},
                { 'type': 'ineq', 'fun': lambda x: x[5] - x22min},
                { 'type': 'ineq', 'fun': lambda x: -x[5] + x22max},
                { 'type': 'ineq', 'fun': lambda x: x[6] - x3min},
                { 'type': 'ineq', 'fun': lambda x: -x[6] + x3max},
                { 'type': 'ineq', 'fun': lambda x: x[7] - x3min},
                { 'type': 'ineq', 'fun': lambda x: -x[7] + x3max},
                { 'type': 'ineq', 'fun': lambda x: x[8] - x3min},
                { 'type': 'ineq', 'fun': lambda x: -x[8] + x3max},
    
    
                )
    
        return cons
    
    
    if __name__ == "__main__":
        # 定义常量值
        args = (41,0.456,0.496,0.524,0.564,47.2,0.456,0.508,58.6,0.458,0.502,0.548)  # Z0,a0,a1,a2,a3,,Z1,b0,b1,Z2,c0,c1,c2
        # 设置参数范围/约束条件
        args1 = (330, 0,25,0,100,0,30,0,50)  # d=700-100-120-150, x1min, x1max,x21min, x21max,x22min,x22max,x3min,x3max
        cons = con(args1)
        # 设置x初始猜测值
        x0 = np.array((10, 10, 10,10,20,20,20,10,10))
        res = minimize(fun(args), x0, method='SLSQP', constraints=cons)
        print('代价',res.fun)
        print(res.success)
        print('解',[np.around(i) for i in res.x])
        
        print('PG1',res.x[0]+res.x[1]+res.x[2]+res.x[3]+100)
        print('PG2', res.x[4] + res.x[5]  + 120)
        print('PG3', res.x[6] + res.x[7] + res.x[8]+150)

![20210220142945889.png_pic_center][]

![在这里插入图片描述][20210220143338101.jpg_pic_center]  
作者：电气-余登武  
![在这里插入图片描述][20210220143424159.gif_pic_center]

[Link 1]: https://blog.csdn.net/KOBEYU652453/article/details/117459596
[2021022009212854.png_pic_center]: /images/20221024/48b62613f43b49a2af0dc07d53f21dc3.png
[20210220092723638.png_pic_center]: https://img-blog.csdnimg.cn/20210220092723638.png#pic_center
[20210220092802271.png_pic_center]: https://img-blog.csdnimg.cn/20210220092802271.png#pic_center
[20210220093528729.png_pic_center]: https://img-blog.csdnimg.cn/20210220093528729.png#pic_center
[20210220093614722.png_pic_center]: /images/20221024/4768a48b77254a74ae586172c71b3c91.png
[20210220093804128.png_pic_center]: https://img-blog.csdnimg.cn/20210220093804128.png#pic_center
[20210220093831278.png_pic_center]: /images/20221024/7b0beb46ea72410f845ad321d9d01766.png
[20210220100525226.png_pic_center]: /images/20221024/c47365dde5104fcdb6b11be91e990035.png
[20210220105800208.png_pic_center]: https://img-blog.csdnimg.cn/20210220105800208.png#pic_center
[2021022010584397.png_pic_center]: /images/20221024/29d2c0f71fbd4a8bbeba3da1888979e2.png
[20210220110600690.png_pic_center]: https://img-blog.csdnimg.cn/20210220110600690.png#pic_center
[20210220110632176.png_pic_center]: /images/20221024/a128a0a8978f46a881a3882ae3d530c7.png
[20210220113624421.png_pic_center]: https://img-blog.csdnimg.cn/20210220113624421.png#pic_center
[20210220113706280.png_pic_center]: /images/20221024/9a97929d187443afbed0ce5aa0375fb7.png
[20210220142945889.png_pic_center]: /images/20221024/8ea74f6097f64c85b2c324c20bb6b95d.png
[20210220143338101.jpg_pic_center]: /images/20221024/7bc0c5b74d054e548cc932f4dcaa81b9.png
[20210220143424159.gif_pic_center]: /images/20221024/b560f687a0304662915123f825535ec6.png