约束二次规划问题

2024-08-08

约束二次规划问题（共5篇）

约束二次规划问题篇1

在制造业中, 制造商追求低成本高效益, 随着制造商掌握越来越多的生产销售信息, 制造商已经开始设计和改进生产运输计划。在这种环境下, 研究有效的生产运输方法得到了学者们的广泛重视。

运输计划通常和选址问题结合起来, 确定工厂和仓库的位置, 使到客户的费用最小化。在选址问题中, 假设货物由工厂运送到仓库, 生产产品的总量和运输的总量是相等的, 用0-1变量表示仓库是否开放, 建立0-1混合整数规划模型。选址问题可以用分支定界法和拉格朗日松弛法解决。

本文研究在生产和仓储过程中的生产运输问题。假设存在多个产销中心, 即生产地同时也是销售地。如果产销中心的需求大于生产能力时, 每个产销中心可以接受来自其它产销中心的产品。通过调节产销中心的产量和运输量, 使得生产和运输中的费用最小化, 制定最优的生产运输方案。

1 问题描述和模型

假设每个产销中心周围存在对产品的需求;产销中心生产每种产品的能力、费用, 及对该种产品确定的需求;在两个产销中心之间有着固定的运输费用;运输所有产品的费用是相同的, 卡车的平均运输能力为M。

整数规划模型如下

$\min \sum_{k = 1}^{n} \sum_{i = 1}^{m} c_{k i} x_{k i} + \sum_{i = 1}^{m} \sum_{\begin{array}{l} j = 1 \\ j \neq 1 \end{array}}^{m} Μ d_{i j} u_{i j} . (1)$

模型参数为:i, j=1, …m为第i、j个产销中心;k=1, …, n为第k种产品;M为运输能力;ski为产销中心i对产品k的需求量;pki为产销中心i对产品k的生产能力;cki为产销中心i生产单位数量产品k所需要的费用;dij为产销中心i运输到产销中心j单位数量产品所需的运输费用。

决策变量定义为:xki为产品k在产销中心i的产量;ykij为从产销中心i运输到产销中心j产品k的数量;uij为从产销中心i运输到产销中心j卡车的数量。

约束条件为

$\begin{array}{l} {\begin{cases} x_{k i} \leq p_{k i}, k = 1, \dots, n, i = 1, \dots, m, \\ x_{k i} \geq \sum_{\begin{array}{l} j = 1 \\ j \neq 1 \end{array}}^{m} y_{k i j}, k = 1, \dots, n, i, j = 1, \dots, m, \\ \sum_{k = 1}^{n} y_{k i j} \leq Μ u_{i j}, i, j = 1, \dots, m, (2) \\ x_{k i} - \sum_{\begin{array}{l} j = 1 \\ j \neq 1 \end{array}}^{m} y_{k i j} + \sum_{\begin{array}{l} j = 1 \\ j \neq 1 \end{array}}^{m} y_{k j i} \geq s_{k i}, k = 1, \dots, n, \\ i, j = 1, \dots, m, \end{cases} \\ {\begin{cases} Μ (u_{i j} - 1) + 1 \leq \sum_{k = 1}^{n} y_{k i j}, i, j = 1, \dots, m, \\ x_{k i} \in Ζ^{+}, y_{k i j} \in Ζ^{+}, k = 1, \dots, n, i, j = 1, \dots, m, \\ u_{i j} \in Ζ^{+}, i, j = 1, \dots, m . \end{cases} \end{array}$

在实际生产过程中, 由于不确定的因素影响, 参数都在一定的区间内变化。假设产销中心i对产品k生产能力为pki, 实际的生产能力由于生产条件和其他客观因素的限制是达不到pki, 实际的生产能力可能是估计生产能力的100 (1-p) %;每个产销中心的需求量也不确定, 它会随着时间的不同及经济形势的变化而围绕着一个值上下波动, 假设产品在每一个产销点的需求量存在一个±s的波动 (100s%) 。鉴于此, 本文将在建立整数规划的同时, 考虑这些不确定因素, 建立一个随机规划模型。

2 随机规划模型的建立

在随机规划模型中引入两个随机参数εpki、εski, 分别代表i产销中心对产品k的实际生产能力和实际需求量, 在一个置信度α的情况下使总费用最小[5,6]。

生产能力和市场需求的约束条件如下:

目标函数为

$\begin{array}{l} z = \sum_{k = 1}^{n} \sum_{k = 1}^{m} c_{k i} x_{k i} + \sum_{i = 1}^{m} \sum_{\begin{array}{l} j = 1 \\ j \neq i \end{array}}^{m} Μ d_{i j} u_{i j}, (3) \\ x_{k i} \leq ε_{p_{k i}}, k = 1, \dots, n, i, j = 1, \dots, m, (4) \\ x_{k i} - \sum_{\begin{array}{l} j = 1 \\ j \neq 1 \end{array}}^{m} y_{k i j} + \sum_{\begin{array}{l} j = 1 \\ j \neq i \end{array}}^{m} y_{k j i} \geq ε_{s_{k i}}, k = 1, \dots, n, i, j = 1, \dots, m . (5) \end{array}$

模型如下:

${\begin{cases} \min \bar{z} \\ Ρ_{r} {ω \in Ω | z = \sum_{k = 1}^{n} \sum_{i = 1}^{m} c_{k i} x_{k i} + \sum_{i = 1}^{m} \sum_{\begin{array}{l} j = 1 \\ j \neq 1 \end{array}}^{m} Μ d_{i j} u_{i j} \leq \bar{z}} \geq α, \\ Ρ_{r} {x_{k i} \leq ε_{p_{k i}}} \geq β_{1}, k = 1, \dots, n, i, j = 1, \dots, m, \\ x_{k i} \geq \sum_{\begin{array}{l} j = 1 \\ j \neq 1 \end{array}}^{m} y_{k i j}, k = 1, \dots, n, i, j = 1, \dots, m, \\ \sum_{k = 1}^{n} y_{k i j} \leq Μ u_{i j}, i, j = 1, \dots, m, (6) \\ Ρ_{r} {x_{k i} - \sum_{\begin{array}{l} k = 1 \\ j \neq 1 \end{array}}^{m} y_{k j i} \geq ε_{S_{k i}}} \geq β_{2}, k = 1, \dots, n, i, j = 1, \dots, m, \\ Μ (u_{i j} - 1) + 1 \leq \sum_{k = 1}^{n} y_{k i j}, i, j = 1, \dots, m, \\ x_{k i} \in Ζ^{+}, y_{k i j} \in Ζ^{+}, k = 1, \dots, n, i, j = 1, \dots, m, \\ u_{i j} \in Ζ^{+}, i, j = 1, \dots, m . \end{cases}$

其中:α、β1、β2为决策者预先给定的置信水平, min $\bar{z}$ 是α的乐观费用。

3 应用算例

5个生产销售中心生产3种不同类型产品a1、a2、a3, 假设运输过程中3种产品可以同时用卡车运输。案例中的生产商拥有5个区域产销中心:Hiroshima、Osaka、Nagoya、Tokyo、Sendai。[7]产品在这5个区域进行运输, 每种产品可以按照形状、颜色、材质分为10到30个类型。本文是为了在它们之间找出生产运输中的最优或满意解。不能取得客户的账单, 因此处理的问题只有5个产销中心和3种产品, 但是这些限制并不影响对问题的研究。每个地方的生产费用、生产能力及需求见表1。

把表1中的参数与随机规划相结合, 假设每个地方市场需求存在20%的上下波动, 每个产销中心的生产能力可能在预计能力的80%, 即s=0.2、p=0.2, εpki, εski 的分布函数为εpki～U (pki, 80%pki) , εski～N (ski, 0.04s $_{k i}^{2}$ ) ) , 置信度α=0.90, β1=β2=0.95。模型的结果见表2、表3。

4 结束语

本文在多个产销中心的前提下研究生产运输规划, 首先运用整数规划建立模型, 然后考虑不确定因素, 并结合算例和随机规划对问题建立随机混合整数规划模型。考虑到参数的随机估计, 最后把模型和生产商的实例相结合, 进行计算机编程求解, 对模型进行测试, 获得满意解。

摘要：从决策者的利益出发对生产销售中心的生产费用和生产销售中心之间的运输费用进行研究使得客户的成本最小化。依据生产销售中心实际生产能力和市场需求, 制定生产销售中心的实际产量和运输量计划, 以达到生产运输成本最小。在不确定环境下, 把随机机会约束和混合整数规划相结合, 计算出在一定置信水平下的最优解, 制定稳健的生产运输方案。

关键词：最优化,生产运输问题,不确定规划,随机机会约束

参考文献

[1]刘宝碇, 赵瑞清, 王纲.不确定规划及应用[M].北京:清华大学出版社, 2003.

[2]刘宝碇, 赵瑞清.随机规划与模糊规划[M].北京:清华大学出版社, 2003.

[3]赵玮, 王荫清.随机运筹学[M].北京:高等教育出版社, 1993.

[4]Liu B D, Iwamura K.Chance constrained programming withfuzzy parameters[J].Fuzzy Set and systems (S1000-1506) , 1998, 94 (2) :227-237.

[5]高雷阜.不确定规划模型及其算法研究[J].辽宁工程技术大学学报, 2003, 22 (3) :413-415.

[6]Masatoshi Sakawa, Ichiro Nishzaki, Yoshio Uemura.Fuzzy pro-gramming and profit and cost allocation for a production and transportation problem[J].EUROPEANJOURNAL OF OPER-ATIONAL RESEARCH, 2000.

约束二次规划问题篇2

产品的拆卸序列规划 (disassembly sequence planning, DSP) 在其数学本质上是一个非确定性多项式 (non-deterministic polynomial hard, NP-hard) 求解的问题。一般而言, 随着产品中零件数量的增加, 产品拆卸方案的求解空间会迅速变得复杂起来, 从而使得采用人工判断的方法难以对其进行有效处理。近年来, 随着人工智能的发展, 越来越多的学者应用启发式算法来解决这一问题, 如Petri网法[1,2]、蚁群算法[3,4]、遗传算法[5,6]以及粒子群算法[7]等。然而, 上述人工智能算法存在陷入局部极值点的问题, 使规划结果不一定为近优解或最优解。约束满足问题 (constraint satisfaction problems, CSP) 作为人工智能领域中一个重要分支, 不存在陷入局部极值点的缺陷, 能较快地找到所有可行解, 通过对可行解进行评价, 即可得到最优的解, 故本文将CSP引入拆卸序列规划中, 通过混合图的指导给出了一种基于CSP的产品拆卸序列规划方法。

1 拆卸序列规划与CSP

拆卸序列规划是指根据产品结构和装配关系等信息, 生成满足一定约束条件的零件拆卸序列的过程[7]。

CSP是一个三元组P={X, D, C}:给定一组有限个变量X={x1, x2, …, xn}, 通过函数变换将每个变量xi一一映射到其对应的有限值域D (xi) ={d1, d2, …, dn}, C为一个有限的约束集合, CSP的解就是在满足所有约束的条件下, 从变量的值域中给每个变量赋值, 使得所有的约束得以满足[8,9,10]。

拆卸序列规划可转化为一类CSP, 拆卸序列规划每一步拆卸的零件可以映射为CSP中的变量X, xi对应拆卸序列中每一步拆卸下来的零件, 拆卸序列规划中每步待拆卸零件可以映射为变量X对应的值域, 零件间拆卸约束关系可映射为CSP中的约束集合C。拆卸序列规划得到的可行序列对应CSP的解。图1给出了产品的拆卸序列规划和CSP的对应关系。

2 基于CSP的拆卸序列规划

2.1 基于CSP的拆卸序列建模

2.1.1 CSP的变量与值域

CSP的解即一个有效的拆卸序列可表示为X={x1, x2, …, xn}, xi为拆卸序列中第i步拆卸的零件。xi的值域di是一个动态变化的域。设A为所有待拆的零件集合, Bi为从第1步到第i步拆卸过程中已拆卸的零件集合, 则第i+1步待拆卸零件xi+1的值域di+1可用下式表示:

$d_{i + 1} = A - B_{i} = \bar{B_{i}} B_{i} \subset A$ (1)

式中, $\bar{B_{i}}$ 为Bi的补集。

2.1.2 CSP的约束

记录产品拆卸信息常用的方法之一是拆卸约束图 (disassembly constraint graph, DCG) [11]。DCG是一个混合图, 可表示为G={V, E, DE} (图2) 。

其中, G表示拆卸约束图;V为节点, 即零件;E为无向边, 表示两个零件间存在连接关系的约束;DE为有向边, 表示两个零件间存在空间干涉, 箭尾零件在空间上阻挡了箭头所指零件的拆卸, 零件的拆卸就是在图中将其对应的节点及与其相关的边删除掉。G可分解为零件间的连接约束图 (混合图中无向边与节点组成的无向图) GE={V, E}和空间干涉图 (混合图中有向边与节点组成的有向图) GDE={V, DE}。设G中节点个数为n, 则图GE和图GDE的邻接矩阵NE和NDE分别为

$Ν_{E} = [\begin{matrix} a_{1, 1}^{(E)} & a_{1, 2}^{(E)} & \dots & a_{1, n}^{(E)} \\ a_{2, 1}^{(E)} & a_{2, 2}^{(E)} & \dots & a_{2, n}^{(E)} \\ ⋮ & ⋮ & ⋮ \\ a_{n, 1}^{(E)} & a_{n, 2}^{(E)} & \dots & a_{n, n}^{(E)} \end{matrix}]$

$\begin{array}{l} Ν_{D E} = [\begin{matrix} a_{1, 1}^{(D E)} & a_{1, 2}^{(D E)} & \dots & a_{1, n}^{(D E)} \\ a_{2, 1}^{(D E)} & a_{2, 2}^{(D E)} & \dots & a_{2, n}^{(D E)} \\ ⋮ & ⋮ & ⋮ \\ a_{n, 1}^{(D E)} & a_{n, 2}^{(D E)} & \dots & a_{n, n}^{(D E)} \end{matrix}] \\ a_{i, j}^{(E)} = {\begin{cases} 0 i 和 j 之间无连接 \\ 1 i 和 j 之间有连接 \end{cases} \end{array}$

且a (E) i, j=a (E) j, i, 当i=j时, a (E) i, j=a (E) j, i=0。

$a_{i, j}^{(D E)} = {\begin{cases} 0 j 在空间上不阻挡 i 的拆卸 \\ 1 j 在空间上阻挡 i 的拆卸 \end{cases}$

且当i=j时, a (DE) i, j=a (DE) j, i=0。

拆卸是一个互相作用的过程, 拆离A、B两个零件, 既可以认为是将零件A从零件B上拆除, 也可以认为是将零件B从A零件上拆除, 故在拆卸序列规划时, 首先需要确定基准件, 认为所有零件均从基准件或基准件所在装配体中拆除, 即拆卸到最后一步产品所剩的零件。同一个装配体在不同的基准件下生成的拆卸混合图极有可能是不同的, 如图3所示。基准件的选取可参照以下原则进行:①在装配体中与其余零件连接关系最多;②体积大, 易保证在拆卸过程中装配体的稳定性;③基准件是拆卸到最后一步所剩的零件, 故其在空间上不能干涉其他零件的拆卸;④方便在拆卸过程中工具和夹具的使用。如在图3所示的装配体中, 以3作为基准件是较为合适的。

一个合理的拆卸序列应保证在操作过程中每一步的零件拆卸不受其他零件干涉, 且零件拆卸后不会影响零件拆卸后剩下的子装配体结构的稳定, 导致拆卸的不确定性, 故在基于CSP的拆卸序列规划中, 约束C可分为两类:连接约束CE和空间干涉约束CDE。

(1) 连接约束CE。

零件拆卸后如剩余子装配体中所有零件均由一定的连接关系组成一个整体, 则认为子装配体是稳定的, 即满足连接约束CE。连接约束CE主要针对图GE进行判断, 当GE去除节点及与其相关的无向边后, 剩余节点及无向边所组成的子图应该是一个连通图, 连接约束CE的逻辑表达如图4所示。

CE的数学表达如下:

if M $_{i i}^{(E)}$ 存在全为零的行

CE不满足

else

CE满足

其中, i为拆卸的零件标识;M $_{i i}^{(E)}$ 为图GE的邻接矩阵NE中元素a (E) i, i (拆卸的零件) 的余子式。若

$Ν_{E} = [\begin{matrix} a_{1, 1}^{(E)} & \dots & a_{1, i - 1}^{(E)} & a_{1, i}^{(E)} & a_{1, i + 1}^{(E)} & \dots & a_{1, n}^{(E)} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ a_{i - 1, 1}^{(E)} & \dots & a_{i - 1, i - 1}^{(E)} & a_{i - 1, i}^{(E)} & a_{i - 1, i + 1}^{(E)} & \dots & a_{i - 1, n}^{(E)} \\ a_{i, 1}^{(E)} & \dots & a_{i, i - 1}^{(E)} & a_{i, i}^{(E)} & a_{i, i + 1}^{(E)} & \dots & a_{i, n}^{(E)} \\ a_{i + 1, 1}^{(E)} & \dots & a_{i + 1, i - 1}^{(E)} & a_{i + 1, i}^{(E)} & a_{i + 1, i + 1}^{(E)} & \dots & a_{i + 1, n}^{(E)} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ a_{n, 1}^{(E)} & \dots & a_{n, i - 1}^{(E)} & a_{n, i}^{(E)} & a_{n, i + 1}^{(E)} & \dots & a_{n, n}^{(E)} \end{matrix}]$

则

$Μ_{i i}^{(E)} = [\begin{matrix} a_{1, 1}^{(E)} & \dots & a_{1, i - 1}^{(E)} & a_{1, i + 1}^{(E)} & \dots & a_{1, n}^{(E)} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ a_{i - 1, 1}^{(E)} & \dots & a_{i - 1, i - 1}^{(E)} & a_{i - 1, i + 1}^{(E)} & \dots & a_{i - 1, n}^{(E)} \\ a_{i + 1, 1}^{(E)} & \dots & a_{i + 1, i - 1}^{(E)} & a_{i + 1, i + 1}^{(E)} & \dots & a_{i + 1, n}^{(E)} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ a_{n, 1}^{(E)} & \dots & a_{n, i - 1}^{(E)} & a_{n, i + 1}^{(E)} & \dots & a_{n, n}^{(E)} \end{matrix}]$

(2) 空间干涉约束CDE。

零件拆卸时, 只有在空间上没有其他零件的阻拦, 才能进行零件的拆卸, 即满足空间干涉约束CDE。空间干涉约束CDE主要针对图GDE进行判断。为满足空间干涉约束CDE, GDE中指向待拆卸零件所表示的节点xi的箭头个数应等于0, 即入度ID (xi) =0。空间干涉约束CDE的逻辑表达如图5所示。

CDE的数学表达如下:

if ID (xi) =0

CDE满足

else

CDE不满足

当某一零件拆卸完成后, 它对子装配体的拆卸约束也随之消失, 因此, 基于CSP的产品拆卸序列规划中约束C是一个动态变化的过程。设xi为第i步拆卸的零件, 则第i+1步待拆零件xi+1的连接约束C $_{E}^{(i + 1)}$ 针对的图G $_{E}^{(i + 1)}$ 的邻接矩阵N $_{E}^{(i + 1)}$ 等于N (i) E中元素a (i) i, i (拆卸的零件) 的余子式M $_{i i}^{(E)}$ ;空间干涉约束C $_{D E}^{(i + 1)}$ 针对的图G $_{D E}^{(i + 1)}$ 的邻接矩阵N $_{D E}^{(i + 1)}$ 等于N $_{D E}^{(i)}$ 中元素a (DE) i, i的余子式M $_{i i}^{(D E)}$ 。

2.2 基于CSP的拆卸序列求解

在确定基准件后, 根据文献[12,13]介绍的方法, 即可通过CAD系统生成拆卸混合图及其相关邻接矩阵, 之后便可采用CSP相关算法进行求解。在众多CSP求解算法中, 回溯算法 (backtracking, BT) 为最常用的算法。该算法通过给部分变量进行复制赋值, 并在赋值过程中对部分赋值的变量进行约束检测, 提前发现不满足约束的变量集合, 显著地提高了搜索效率[14,15], 故本文采用回溯算法针对拆卸序列规划进行求解。由拆卸的顺序性可知, 当拆卸到装配体只剩两个零件A、B (B为基准件) 时, 其序列只能为A→B。因此, 为减小运算量, 在求解拆卸序列时, 首先对基准件赋值, 再使用回溯算法对拆卸的第一个零件到倒数第三个零件的序列进行求解, 拆卸序列中倒数第二个零件则在回溯算法求解后自动得到。当得到可行拆卸序列后, 设计人员即可根据具体要求对可行拆卸序列进行评价, 得到最符合要求的拆卸序列。由于评价得到的结果是从所有可行解中得到的, 容易证明, 得到的拆卸序列应是全局最优解。以拆卸时间最小为例, 则对拆卸序列的评价函数为

$\min X = \sum_{i = 1}^{n} Τ_{Ρ} (x_{i}) + \sum_{i = 1}^{n} Τ_{D} (x_{i}) + \sum_{i = 1}^{m} Τ_{C} (i) (2)$

式中, TP (xi) 为拆卸零件xi之前的准备时间, 主要包括更换拆卸工具消耗的时间以及将拆卸工具放置到拆卸位置的时间;TC (i) 为拆卸零件过程中每次对装配体旋转、移动以便进行拆卸所消耗的时间;m为拆卸过程中对装配体旋转、移动的次数;TD (xi) 为拆卸零件xi的工作时间。

基于CSP的拆卸序列求解的流程如图6所示。

基于回溯算法的产品拆卸序列求解是一个递归算法, 即

3 应用实例

以一个简化的洗碗机门体作为应用实例, 如图7所示 (为减小运算量, 部分细小零件以及螺钉等连接件已经移除) 。产品采用手工拆卸的方式, 目标为拆卸时间最短。根据基准件选取原则, 采用“4-内门板”为基准件, 生成的洗碗机门体拆卸混合图见图8。洗碗机门体零件拆卸信息如表1所示。由图8得到洗碗机门体连接约束图GE的邻接矩阵为

$Ν_{E} = [\begin{matrix} 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 & 1 & 1 & 1 \\ 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \end{matrix}]$

空间干涉图GDE的邻接矩阵为

$Ν_{D E} = [\begin{matrix} 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \end{matrix}]$

1.控制面板组件 2.锁扣夹板 3.门锁扣 4.内门板 5.分配器组件 6.旋钮 7.门铰链组件 8.门保护条 9.外门板

注:零件1、9虽然拆卸方向为+X, 但由于其螺钉连接固定在-X方向, 故其拆卸应先从-X将固定的螺钉去除, 然后再从+X取走零件, 零件5的情况与零件1、9正好相反。

选择一次运算情况对算法进行说明。首先对基准件进行定义, 确定x9=4。根据本文所述算法进行第一次递归调用, 运行函数BT (1) , 令x1=1, 对C (1) E以及C (1) DE进行判断, 确定满足约束。之后进行第二次递归调用, 运行函数BT (2) , 令x2=2, 对C (2) E以及C (2) DE进行判断, 确定满足约束。如此进行至第四次递归, 运行函数BT (5) (此时x2=2, x3=3) , 令x4=5, 发现同时不满足约束C (4) E和C (4) DE, 则不再对剩余变量进行赋值, 开始回溯, 对x4重新赋值为6。发现此时满足约束, 进行第5步递归调用。最后一直运行到第7步递归调用, 得到一个解{x1, x2, x3, x4, x5, x6, x7}={1, 2, 3, 6, 7, 9, 5}。令x8为最后与基准件相连零件8, 最终得到一个可行拆卸序列X={x1, x2, x3, x4, x5, x6, x7, x8, x9}={1, 2, 3, 6, 7, 9, 5, 8, 4}。编程对问题进行求解, 即可得到洗碗机门体所有可行拆卸序列。

当确定所有可行拆卸序列后, 即可对可行拆卸序列进行评价, 得到最符合要求的拆卸序列。本例采用手工拆卸的方式对洗碗机门体进行拆卸, 拆卸过程中洗碗机门体平放在工作台上, 初始状况默认图7中+X方向向上。在拆卸过程中每次进行翻转的时间TC (i) 近似相同约为8s。由于洗碗机门体较小, 因此将拆卸工具放置到拆卸位置对拆卸影响也较小, 平均约为1s, 而更换拆卸工具的时间为4s, 故TP (xi) =1+4=5s。根据式 (2) , 对得到的所有可行解进行评价, 结果表明, 洗碗机门体所有可行拆卸序列在上述条件下拆卸耗时均近似相同, 表2给出了部分拆卸时间最小的序列。

4 结语

本文在拆卸混合图的基础上, 将拆卸序列规划转化为一类约束满足问题 (CSP) 进行求解, 给出了基于CSP的拆卸序列求解模型及相关算法, 并用实例进行了说明。相较基于蚁群算法、遗传算法以及粒子群算法等人工智能算法的拆卸序列规划而言, 基于CSP的拆卸序列规划不存在结果陷入局部极值点的缺陷。同时, 基于CSP的拆卸序列规划不对生成的序列进行评价、优选, 而是在得出所有可行序列后另行评价优选, 故当评价指标改变后, 只需对CSP的解重新进行评价优选, 而不必再次进行整体的搜索, 因而提高了效率和适应性。

约束二次规划问题篇3

随着全球能源和环境问题的日益突出,风能等可再生能源发电技术得到迅速发展,风电并网的规模也越来越大[1,2]。由于风电出力具有很强的不确定性,含风电场的电网日前发电调度问题常描述成为一个含有随机变量的动态经济调度(DED)问题[3,4]。为了使获得的发电调度计划对于风电场出力不确定性具有适应性,通常采用场景法,通过对风电场出力随机变量进行抽样模拟,进而将随机模型转换为确定性DED模型[5,6,7,8,9,10]。由于随着抽样的场景数目的增多,场景法求解随机DED问题的模型维数将快速增大,直接求解非常费时,因而目前该方法主要应用于中小型系统的优化调度,应用于实际大型电网将面临模型维数太大、求解时间太长的问题。

另一方面,由于发电机组相邻时段出力的变化量存在爬坡率的限制,含风电场的电网日前发电调度问题是一个含有一天所有时段变量的联合优化模型,所有时段变量的同时求解是导致问题维数太大的另一个关键因素。动态规划(DP)法根据最优性原理,即Bellman方程可实现对于日前发电优化调度问题各个时段决策量的解耦求解[11,12]。然而,实际大电网机组众多,每个时段各个机组出力组成的状态维数非常之大,DP法应用于大电网发电调度问题将不可避免地面临着“维数灾”难题。

近似动态规划(ADP)理论通过近似描述值函数与状态量之间的关系来克服“维数灾”难题,文献[13,14]应用ADP理论求解大规模机组组合(UC)问题,不过没有考虑风电随机性对于电网UC的影响。文献[15,16]将ADP理论应用于含风电和储能装置的小型系统优化调度。文献[17]将含有单一风电场和抽水蓄能电站的电力系统随机DED问题描述为随机型存储模型,以抽水蓄能电站水库的储水电量作为系统的存储水平,并采用ADP算法克服随机规划问题中目标函数含有数学期望计算的难题。然而,所提方法只适用于必须含有抽水蓄能电站的电网调度问题;且建立的模型中并没有考虑网络安全约束,获得的调度计划无法满足工程应用需求;另外,目标函数采用机组出力的一次函数,能否适应于DED问题通常采用的二次目标函数还有待进一步验证。

由于目前国内大部分省级电网中不含有抽水蓄能电站,对于不含有抽蓄电站的大型电网,如何应用ADP算法求解其随机DED问题,同时考虑网络安全约束的影响,对于扩大ADP算法在求解随机优化调度问题方面的适用范围,无疑具有重要的实用意义。因而,本文以系统中多个风电场出力的日前预测曲线为基础场景,借助拉丁超立方抽样生成风电场出力误差场景。以当前时段的系统正旋转备用容量作为资源存储量,列写了相邻时段关系的系统状态转移方程,从而建立了不含抽水蓄能电站电网的随机DED问题的随机型虚拟存储器模型(VSM)。在考虑网络安全约束的条件下利用误差场景对随机DED问题各个时段的值函数进行训练,利用训练得到的值函数对预测场景下的VSM进行求解,得到考虑风电出力随机性影响的常规机组日前发电出力计划。

1 随机型VSM描述

存储模型通过设置一个表示系统资源存储量的变量作为系统的状态变量,很好地解决动态规划问题状态的“维数灾”。由文献[17]可知,对于含有抽水蓄能电站的电网,可以方便地以抽水蓄能电站水库的储水电量作为系统的资源存储量,但对不含抽水蓄能电站的电网,在系统中难以找到可直接表示系统资源存储量的变量,因此如何选取系统的资源存储量,成为此类系统存储模型构建的难点和应用ADP算法求解该类系统随机DED问题的关键。

由于在一般电力系统中,系统的旋转备用容量反映了系统的可调控发电能力,相当于存储在系统中可用于平衡风电场出力随机波动和负荷需求变化的容量,由于存储模型只设置一个表示资源存储量的变量,故本文选取系统的正旋转备用容量作为系统的资源存储量,并根据相邻时段的系统正旋转备用容量的变化关系,列写出系统的状态转移方程,从而建立适用于一般电力系统随机DED问题的VSM,并采用ADP算法求解。

1.1 目标函数

优化目标取常规机组总发电燃料耗量最小,由于风电出力的随机性,目标函数应表示为风电的各种可能出力下对应的常规机组总发电燃料耗量的期望值最小,如式(1)所示。

式中:T为调度周期总的时段数,本文取T=96;ΔT为每个时段的持续时间,即15min;St为t时段系统所处状态;xt为决策变量向量;Ct(St,xt)为时段t所有NgNg台常规机组的燃料耗量,

,其中,Pi,t为第i台常规机组在时段t的发电功率,Ai,2,Ai,1,Ai,0为第i台常规机组的耗量特性系数,对于水电机组,有Ai,2=Ai,1=Ai,0=0;E{·}为期望函数;Πt为xt的可行域。

1.2 约束条件

1)基本约束

式中:Ploadj,t为负荷节点j在时段t的功率预测值;Nd为负荷节点数;Pwk,t为风电场k在时段t的出力值,为随机变量;Pi-和P-i分别为机组i的有功出力上、下限;rdi和rui分别为机组i的向下、上爬坡率。

其中,第1个式子为功率平衡方程,第2个式子为常规机组的有功出力上、下限约束,第3个式子为常规机组的爬坡约束。

2)网络安全约束

式中:Fl,t为时段t第l条支路的传输功率;Flmin和Flmax分别为第l条支路传输功率的下限和上限,一般Flmin直接取Flmax的负值;Fij,t为第i个安全断面中第j条支路在时段t的传输功率;Ωi为第i个断面包含的支路集合;FΩimin和FΩimax分别为第i个断面的安全下限和上限。

其中,第1个式子为输电线路安全约束,第2个式子为断面安全约束。支路传输功率Fl,t可由直流潮流模型近似表示为:

式中:Gl,i,Dl,j,Wl,k分别为第i台常规机组、第j个负荷和第k个风电场对支路l的功率转移分布因子,其值由网络结构和支路参数确定[18]。

由于实际大电网支路数众多,若在模型中加入所有的支路安全约束,优化模型的规模会大幅度增加,进而导致求解速度的快速下降。本文采用“求解→安全校验→添加越限支路约束再求解”循环的方法,直至所有支路功率都通过安全校验,这样可加快求解速度,并得到满足所有支路安全约束的最优解[19]。

3)旋转备用约束

为应对风电出力的不确定性和负荷预测误差,系统中应保留一定的旋转备用容量以保证系统安全可靠运行。系统及各常规机组备用约束如下:

式中:sui,t和sdi,t分别为机组i在时段t能够提供的正旋转备用容量和负旋转备用容量;T10为要求的机组旋转备用响应时间,取10min;Su,t和Sd,t分别为系统在时段t的正、负旋转备用容量;Lu和Ld分别为负荷对系统正、负旋转备用的需求系数,通常设定为2%~5%;wu和wd分别为风电场出力对系统正、负旋转备用的需求系数,根据目前国内风电功率预测系统的预测误差范围,可设定为10%~25%;P-wk为风电场k的额定出力。

4)状态转移方程约束

通过将系统正旋转备用容量Su,t设置为系统在时段t的资源存储量,取系统时段t的状态向量为St=(Su,t,Pw,t),则系统的状态转移方程如下:

式中:Ps,t为时段t系统正旋转备用容量相对上一时段的变化量;Pw,t为时段t所有风电场出力组成的向量。Ps,t既与时段t风电场出力随机变量Pw,t有关,又与时段t常规机组出力决策变量xt有关。该方程的物理意义是系统状态在随机变量和决策变量共同作用下的演化形式,体现了相邻两个时段系统正旋转备用容量之间的耦合关系。

5)系统旋转备用变化量约束

每一时段系统正旋转备用容量相对上一时段的变化量有一定的范围限制,这个范围可由风电出力变化量与负荷变化量确定。当风电出力增加大于负荷增长时,系统正旋转备用变化量应满足:

当风电出力增加小于负荷增长时,系统正旋转备用变化量应满足:

2 ADP思想与模型处理

2.1 DP的局限性

基于Bellman的最优性原理,求解多阶段决策问题时,严格意义上DP可以求得全局最优解[20]。对初值问题DP的求解决策过程如图1所示。图中:Jt为时段t的收益;St=fs(St-1,xt)为时段t-1到时段t的状态转移方程。令xt*为时段t的最优决策,求解时先从最后时段开始往前逐一时段递推,依次得到各时段最优决策和值函数与状态关系xT*(ST),VT(ST),x*T-1(ST-1),VT-1(ST-1),…,x1*(S1),V1(S1)的表达式,其中,Vt为时段t的值函数,即从时段t到末时段T内所有阶段收益总和的最优值,然后代入初始状态S0并结合状态转移方程,从前往后逐一求得各时段的最优决策和值函数。

由DP的求解过程可以看出,应用DP求解DED问题,当机组出力连续时,由于爬坡率约束的存在,相邻时段之间的决策变量具有耦合,机组出力可行域也是随不同时段变化的,难以用解析表达式描述决策、收益与状态之间的关系;当机组出力离散时,可以对所有的机组出力组合情况进行评估,但随着机组数、时段数、状态变量数的增加,组合情况呈指数式增长,将面临“维数灾”问题。

2.2 ADP思想

由DP的决策过程可知,DP在求解DED问题时虽然能够求得全局最优解,但对于实际大型电网来说其推导过程过于繁琐,求解的复杂程度难以接受。近年来,Powell等人将ADP方法运用到具有随机性可再生能源接入的电力系统调度中,很好地克服了DP求解DED问题的局限性。

由2.1节可知,DP在决策前需从后往前逐一推导每一状态St对应的值函数Vt(St)的表达式,这是DP求解的关键和难点。如果假定各时段值函数的表达式Vt(St)已知,则在求解当前时段t时,只需在St-1的基础上结合状态转移方程St=fs(St-1,xt)和当前时段值函数Vt(St),即可求得当前时段t的最优决策xt*。但各时段值函数的精确表达式Vt(St)事先无法预知,这为模型的解耦求解带来困难,ADP的思想就是通过采用近似值函数来逼近描述时段t的值函数与状态St的关系,从而实现模型的时段解耦求解,进而可依次求得各时段的近似最优决策xt。由此可以看出,ADP算法的关键就是近似值函数的合理表示。

2.3 模型处理

为了方便应用ADP对随机DED问题的VSM进行求解,需对模型进行一些必要的处理。为此将每个时段假想成两个阶段,分别对应决策前状态(Su,t,Pw,t)和决策后状态(Sxu,t,Pw,t)[21],并定义S^u,t(Pw,t)为时段t观察到随机变量的实现值后状态的变化量。其中,决策前状态(Su,t,Pw,t)表示仅考虑随机变量引起的状态变化量S^u,t(Pw,t)的作用,而未做出决策前的系统状态;决策后状态(Sxu,t,Pw,t)表示做出最优决策后系统的状态。因此系统状态转移方程转化为:

式(9)表示假定时段t观察到的风电变化量直接作用于系统正旋转备用容量,由Sxu,t-1增加演化为Su,t;式(10)表示做出决策得到常规机组出力值xt后,Su,t加上系统正旋转备用容量的实际变化量Ps,t(xt),并扣除没有实际作用效果的后,最终得到决策后系统正旋转备用容量Sxu,t。引入决策前状态和决策后状态后,可得时段t的决策前状态值函数Vt*(Su,t,Pw,t)和决策后状态值函数Vtx(Sxu,t,Pw,t)如下:

此处Πt为由式(2)至式(5)和式(7)、式(8)所确定的xt的可行域。

由式(9)可知,从时段t的决策后状态到时段t+1的决策前状态,仅考虑随机因素的作用,所以式(12)中含期望计算,这给求解带来不便。因此在应用ADP算法求解随机DED问题时,除了要解决近似值函数的合理描述问题,还要处理好系统中随机因素引起的期望计算。

根据文献[21]可知,对于资源分配问题,对于没有明显特性的值函数,可以通过查表与聚类、参数模型、非参数模型等一般工具获得近似值函数;而对于值函数相对资源存储量具有连续、线性或近似线性、非线性(凹性或凸性)性质的,可以采用接近其特性的函数对值函数进行近似。文献[22]给出了对于线性目标函数存储模型采用满足凸性的分段线性函数近似值函数的收敛证明,由于上述VSM的目标函数是二次函数,和线性函数一样具有凸函数特性,因而本文采用满足凸性的分段线性函数来逼近其决策后状态的值函数Vtx(Sxu,t,Pw,t)。因此,通过在决策后状态Sux,t的取值区间上取离散断点R=ρ,2ρ,…,mρ,令vt(Pw,t)=[vt(Pw,t,ρ),vt(Pw,t,2ρ),…,vt(Pw,t,mρ)]T为时段t值函数的斜率向量,其中,m为存储量的所有分段数,ρ为每段长度,则t时段决策后状态的近似值函数可表示如下:

式中:Vtb为时段t值函数的截距;ytr为第r段的存储量。

将式(13)代入式(11),则随机DED问题的VSM可转化为如下不含期望运算的确定性二次规划模型:

3 VSM的ADP求解

3.1 近似值函数的求取

应用ADP求解VSM时,近似值函数t(Sxu,t,Pw,t)对精确值函数Vtx(Sxu,t,Pw,t)的近似精度越高,则近似最优决策xt越接近xt*。为获得高质量的近似值函数,首先根据确定性优化模型求解结果对近似值函数的斜率向量和截距进行初始化,然后扫描误差场景,在每个场景下逐个时间段求解二次规划问题(式(14)),并根据求解结果采用逐次投影近似路径(SPAR)算法[16]修正每次迭代的近似斜率值vtn(Pw,t)和截距值Vntb,直到得到收敛的近似值函数tn(Sxu,t,Pw,t)。SPAR算法对近似值函数的求取过程如图2所示。图中,tn(Sxu,t,Pw,t)为第n次迭代所得近似值函数,Vtx(Sxu,t,Pw,t)为事先未知的精确值函数,和vtx分别为第n次迭代时段t值函数的斜率近似值和时段t值函数斜率的精确值。

斜率向量和截距初始化时,首先根据确定性优化模型的决策结果,获得各时段的资源存储水平Sux,,t0及相应时段的值函数值Vt0。斜率初值设定时以(Sux,,t0,Vt0)作为该时段值函数的极小值点,且其两边各段的斜率符号相反,与极小值点相邻的两段关键点的斜率初始值可根据优化目标的物理意义合理设定,本文主要根据常规机组的煤耗特性系数确定,其余段的斜率根据满足值函数凸性的斜率单调递增特性依次设定。在初始斜率向量给定后,初始截距V0tb根据式(15)确定。

式中:为时段t值函数的斜率初始值。

给定初始斜率向量和截距后,依次在每个场景下逐个时段求解二次规划模型(式(14)),再进行斜率和截距修正,斜率修正过程参见文献[17],得到第n个场景迭代的近似斜率分量和近似值函数值tn(·,Pnw,t)后,根据式(16)计算截距修正值Vntb为:

实际计算中,可只对图2所示关键区域的两段斜率进行修正,再结合截距修正,以节省值函数训练时间,提高计算速度。

3.2 ADP求解过程

ADP求解随机DED问题VSM的步骤如下。

步骤1:求解预测场景对应的确定性经济调度模型,得到各时段决策xt0、存储量和值函数值Vt0。

步骤2:初始化各时段的近似斜率向量,根据初始斜率值和来确定初始截距V0tb。

步骤3:借助拉丁超立方抽样生成基于预测场景P0w,1,P0w,2,…,P0w,T的误差场景样本,获得N个误差场景Pnw,1,Pnw,2,…,Pnw,T(n=1,2,…,N)[23];令n=1,t=1。

步骤4:若n>N则转步骤11,否则继续。

步骤5:若t>T则转步骤9;若t=1,则令的上限和下限设置为;否则计算决策前的资源存储量

步骤6:求解式(14)的二次规划模型,得到最优决策xtn,并计算得到决策后的资源存储量

步骤7:若t<T,则进行斜率和截距修正。

步骤8:t增加1,转步骤5。

步骤9:对场景n的求解结果进行网络安全校验,若存在支路越限,则将越限支路的安全约束加到式(14)所示模型,令t=1,转步骤5;若不存在支路越限,则转步骤10。

步骤10:n增加1,转步骤4。

步骤11:求解预测场景的VSM,获得调度计划。

4 算例分析

为验证本文所建立的随机DED问题的VSM和ADP求解算法的有效性,对某个不含抽水蓄能电站省级电网的发电调度进行建模和求解。以该省网2015年1月5号的数据为例,共有85台常规机组,其中火电机组46台,装机容量为14 560 MW;水电机组39台,装机容量为8 208 MW。风电场5座,额定容量分别为3 958.5,1 140,192,99,49.6 MW,其并网站点见附录A图A1,其中前两个风电场的出力预测曲线,以及系统日前负荷预测曲线和外送功率曲线见附录A图A2和图A3。系统共有线路498条,3个安全断面,各断面数据见附录A表A1。

假定风电出力预测误差服从正态分布,数学期望为各时刻的风电出力预测值,标准差为预测值的20%,借助拉丁超立方抽样方法分别生成20,50,100,200个误差场景进行求解。以20个场景的求解为例,训练过程中值函数变化如图3所示。可以看出,训练刚开始时误差场景的值函数与由确定性模型优化结果反推的值函数非常接近,随着训练的进行,后面误差场景求解得到的值函数慢慢趋向收敛,整个训练过程耗时198.39s。

本文构建的随机VSM和ADP算法求解结果与场景法求解结果的值函数对比见附录A图A4。采用本文模型和ADP算法求得的一天总发电燃料耗量为7.572 027万t,场景法求得的总发电燃料耗量为7.487 056万t,且由附录A图A4中各时段的值函数比较可以看出,ADP算法与场景法求得的燃料耗量结果十分接近。以上比较充分说明了本文建立的不含抽水蓄能电站的随机DED问题的VSM及ADP算法求解的正确有效性。

ADP算法求得的系统正旋转备用与场景法优化结果比较如图4所示。可以看出,两种方法得到的系统正旋转备用的整体变化趋势也基本一致,只是ADP算法得到的系统正旋转备用整体上比场景法略微大一些。

两种方法得到的机组出力计划比较如图5和图6所示。由图5可以看出,两种方法得到的火电机组的出力计划基本一致,部分机组在某些时段出力存在微小偏差。由图6可以看出,场景法得到的水电机组出力存在很大的跳跃,而ADP算法得到的水电机组出力则变化比较缓慢,这是由于水电机组功率调节速度快,每个时段可调节功率范围较大,因此场景法求解时在满足各种约束的条件下为了优化目标函数而使得机组出力会有较大的波动跳跃,这与水电机组自身的调节特性相吻合,而在采用VSM和ADP算法求解时由于式(6)至式(8)的约束,限制了系统正旋转备用的变化,使得备用响应容量较大的水电机组的出力变化也较为缓慢,这更符合实际电网运行调度中对机组出力的调控要求。

同时,由于模型中添加了断面安全约束,能够保证所获得调度方案下系统的安全运行。以20个误差场景的优化为例,与不含断面安全约束求解结果对应的安全断面2的输电功率对比如表1所示。可以看到,在未加断面约束时优化得到的总燃料耗量为75 706.61t,但断面2在某些时段存在功率越限;加入断面约束后,总燃料耗量为75 720.27t,比不加断面约束时增加了13.66t,但断面2功率都小于安全极限。因此,在模型中加入网络安全约束后,为了使系统的关键线路和断面的输送功率在限定范围内,机组的出力安排可能会使得系统总的燃料耗量有所增加,这在一定程度上使得系统的经济效益有所下降,但却避免了系统运行在不安全状态,对系统的安全可靠运行具有重要意义。

接下来分别将该算法与场景法在20,50,100,200个场景的情况下进行比较,验证该算法的计算性能。使用计算机为Intel(R)Core(TM)i7-4900MQ CPU 2.80GHz/32GB内存,计算结果如表2所示。由表2可见,场景法在场景数较少时具有较快的计算速度,但随着场景数的增加,计算所需内存和时间都大幅增长,这在很大程度上限制了场景法的应用,尤其是对于风电场数目多需要抽样很多个场景来准确模拟风电出力特性的大型电网调度问题,场景法求解将会受到计算机内存容量限制。而ADP算法由于实现了对各个场景和各个时段的解耦求解,将大规模优化问题分解成若干个小规模优化问题逐个求解,所以随着场景数的增加,所需内存无明显增长,求解时间也基本只增加了新增加场景进行值函数训练所增加的时间。对于100个场景求解时间只有16min左右,约为场景法的1/12;即使对于200个场景求解时间也只有33min左右,计算速度明显提高。

同时,将所提出算法与基于极限场景集的鲁棒优化调度(RS)方法比较[24]。为保证极限场景能覆盖95%的可能风电出力,取风电功率的变化范围为[μ-2σ,μ+2σ],其中,μ为期望值,σ为标准差值,由于系统中含有5个风电场,故共有25即32个极限场景,RS方法求解总耗时6 378.83s,优化结果的总燃料耗量为75 654.04t。

由此可以看出,虽然RS方法比场景法更能保证对风电出力大范围波动的适应性,但其目标函数值也更大,且在极限场景只有32个的情况下,其求解时间已经分别达到50个场景下场景法和ADP算法的3.3倍和12.9倍,当系统中风电场数目增大时,其求解时间将增加得更为明显。因此,ADP算法与RS方法比较同样能够大幅提高计算速度。

另外,由于极限场景的数目与风电场数目呈指数关系增长,随着风电场数目的增大,RS方法和场景法一样会面临由于问题规模过大超出计算机内存容量限制进而无法求解的问题。因此,ADP算法对于含多个风电场的大型电网随机优化调度问题具有更好的适应性,在求解速度上相对RS方法及场景法具有明显的优势,能够很好地满足应用于实际大型电网日前发电调度的要求。

5 结语

本文将ADP理论推广应用于不含抽水蓄能电站的电网随机DED问题,以正旋转备用容量为存储量,建立不含抽水蓄能电站的电网安全约束随机DED问题的VSM,并通过与场景法和鲁棒优化调度方法求解结果的比较分析验证了所建模型和求解算法的正确有效性,为ADP理论应用于快速求解一般大型电网的随机DED问题提供了新途径。ADP算法实现了对随机优化调度模型各个场景和各个时段的解耦求解,将一个大规模优化问题分解为一系列小规模优化问题,有效提高了对大电网随机优化调度模型的求解速度。采用ADP算法求解随机型VSM的优化结果中对应的水电机组出力变化比场景法更加合理,符合实际电网运行调度中对机组出力的调控要求。另外,对于含有抽水蓄能电站的电网调度问题,也可以采用本文提出的VSM建模方法并通过ADP算法快速求解;即便是对于含有多个抽水蓄能电站的电网调度问题,文献[17]的建模方法由于只适用于含单一抽水蓄能电站的电网,会存在建模困难,而本文的VSM建模方法同样能够适用。

本文研究中采用分段线性函数对值函数进行近似,所得调度方案对应的目标函数值比场景法的结果有所增大,如何提高值函数的近似精度,以获得更优的调度方案是本文下一步工作重点;同时,本文建立模型中未考虑不同时段机组启停状态的变化,如何应用ADP算法求解随机机组组合问题是本文的进一步研究方向。

附录见本刊网络版(http://www.aeps-info.com/aeps/ch/index.aspx)。

摘要：针对大电网安全约束随机动态经济调度(DED)问题的求解时间太长,提出了应用近似动态规划算法快速求解不含抽水蓄能电站电网的安全约束随机DED问题的方法。建立了随机DED问题的虚拟存储器模型,以系统的正旋转备用容量作为存储变量,构建系统相邻时段的状态转移方程,并考虑了各输电线路和断面的安全约束。以风电场日前功率预测曲线为基础,通过拉丁超立方抽样产生风电场出力的误差场景,并逐一场景递推求解每个时段的二次规划模型以对各个时段的值函数进行训练,形成收敛的值函数,再代入预测场景求解以获得最终的优化调度方案。该方法实现了对随机DED模型各个场景和各个时段的解耦求解,将一个大规模优化问题分解为一系列的小规模优化问题,有效提高了对大电网随机DED模型的求解速度。以某一实际省级电网为算例,通过与场景法和鲁棒优化调度方法的比较验证了所提出模型和求解方法的正确有效性。

约束二次规划问题篇4

通常, 集装箱码头岸桥调度问题 (QCSP:Quay Crane Scheduling Problem) 是研究在一定约束下寻求一个最优的岸桥调度计划以使得整个船舶的装卸时间最短。在岸桥调度问题中, 一组连续的装卸操作称为一个任务, 不同文献对任务粒度的约定不同。任务可以定义在:一个贝区域、一个完整的贝、同一贝中的一组集装箱、或单个集装箱上, 参见Bierwirth和Meisel (2009) [1]。其中以同一贝中的一组集装箱作为一个调度任务 (一个贝中可以有多组任务) 兼具灵活性和便利性, 本文的模型即采用此定义。

Daganzo (1989) [2]较早研究QCSP问题, 他提出了一个MIP模型, 但模型未考虑岸桥冲突。针对该模型, Peterkofsky 和Daganzo (1990) [3]提出了一个改进的分支定界算法进行求解。此后岸桥调度问题的研究得到越来越多的关注。

在以同一贝中的一组集装箱作为一个调度任务的QCSP问题中, Kim和Park (2004) [4]提出了一个考虑任务顺序关系及岸桥冲突约束的MIP模型, 并设计了分支定界算法 (B&B:Branch-and-Bound) 和贪婪随机适应性搜索算法 (GRASP:Greedy Randomized Adaptive Search Procedure) 求解, 该模型是很多后续模型的基础。Moccia等 (2006) [5]发现Kim的模型在可能会产生一些岸桥冲突的解, 提出了一个修正的MIP模型, 并设计了分支割算法 (B&C:Branch-and-Cut) , 改进了解的质量。Sammarra等 (2007) [6]进一步设计了禁忌搜索算法 (TS:Tabu Search) 求解该修正模型, 求解时间大大缩短, 但解的质量稍差。Bierwirth和Meisel (2009) [1]发现Moccia的模型在某些情况下仍然会产生岸桥冲突的解, 为此又提出了一个修正的MIP模型解决岸桥冲突, 并设计了一个基于分支定界的单向调度启发式算法 (UDS:Unidirectional Scheduling) 求解, 算例表明该算法能够在更短时间内产生比以往文献质量更好的解。Chung和Choy (2012) [7]设计了一个改进的遗传算法GA求解QCSP, 其质量和效率未超过Bierwirth和Meisel (2009) [1]。

有些文献不以同一贝中的一组集装箱作为任务调度单元, 这方面的研究主要有Lim等 (2004, 2007) [8,9], Zhu和Lim (2006) [10], Liu等 (2006) [11]。

国内, 曾庆成和高宇 (2006) [12]将贪婪随机适应性搜索算法结合进遗传算法对Kim和Park (2004) [4]的模型进行算例实验求解, 对中小规模的问题给出了求解结果, 未与国外现有结果进行性能比较。董良才等 (2011) [13]提出了带时间窗的岸桥调度MIP模型, 设计了遗传算法, 该模型考虑了岸桥交叉干扰, 但未考虑岸桥最小安全距离约束。李晨等 (2010) [14]设计了一种遗传算法求解岸桥调度问题, 但该文以跨越多贝的区域为任务单元。韩骏等 (2006) [15]研究了集装箱码头泊位与岸桥协调调度优化问题, 其岸桥调度部分对岸桥不能互相跨越仅用同一艘船只能使用连续岸桥表示, 不能够约束所有岸桥跨越情况。

在以同一贝中的一组集装箱作为任务调度单元的算法中, Bierwirth和Meisel (2009) [1]的UDS算法是质量和效率较好的算法, 但该方法都仅能求得岸桥单向移动 (调度) 解, 即其解中所有岸桥都只能向同一个方向移动。本文将设计一个新的约束规划模型, 理论上可以求得完全的岸桥双向移动解, 即任一岸桥都可以双向移动。

本文的主要贡献在于提出了一个能够求得岸桥双向调度解的约束规划 (CP:constraint programming) 模型, 通过数值实验验证其有效性。CP模型中模型表示与求解算法分离, 易于修改约束, 实验结果显示求解效率也较好, 取得了模型灵活性和求解效率较好的平衡。

1 问题描述

集装箱船到达码头后, 要进行集装箱装卸操作。集装箱船上的集装箱是按贝存放的, 可以根据一定规则将同一贝中的相邻集装箱分成若干组, 例如可以按同一目的港、同一箱型 (对装船而言) 分组, 或按同一装货港、同一箱型 (对卸船而言) 分组。一个调度任务即为对一组集装箱的装 (或卸) 载操作, 不能中断。

码头通常分配多个岸桥同时装卸一条船, 岸桥在同一轨道上运行, 因此不能互相跨越, 并且岸桥间必须保持一个安全间距, 这就是岸桥冲突约束, 它是QCSP问题的一个关键约束。同一贝的任务之间可能存在顺序关系约束, 比如卸载任务必须早于装载任务, 甲板下的装载任务必须早于甲板上的装载任务等。此外, 有些任务不能同时操作, 例如多个装载任务需要从堆场同一箱区提箱, 可能导致堆场该区域拥堵, 这些装载任务就不能同时操作。岸桥调度就是遵守这些约束的情况下为各岸桥分配任务以及各任务的开始结束时间, 以使得整船装卸完工时间最短。本问题详细描述参见Bierwirth和Meisel (2009) [1]。

2 QCSP约束规划模型构建

2.1 模型特征

CP模型在不同CP系统中表示方法不同, 本文所建CP模型是用IBM ILOG CPLEX Optimization Studio 12.4中的OPL语言实现的, 使用了OPL的语言构造, 这些构造许多都是含义自明的, 以下必要时会在文中出现时解释其含义。OPL针对调度问题设计了一个特别的数学概念, 即区间变量, 以简化模型表示和提升搜索效率, 本文的模型也是以区间变量为核心设计的, 区间变量的含义在2.2节 (3) 介绍。

2.2 模型符号

(1) 集合/索引

Tasks/i, j : 装卸任务, |Tasks|=n

Q/k, v, w: 岸桥, 岸桥沿着船贝的方向从左到右递增编号, |Q|=q

Phi:具有顺序关系的任务对集 (i, j) 构成的集合, i必须在j前执行, 这里i, j同贝

Psi: 由于某种外部因素导致不能同时执行的任务对 (i, j) 构成的集合

transitionTimes: 自定义集合, 用于表示岸桥在不同贝之间移动的时间, 集合表达式为{ (b1, b2, |b1-b2|·t) |b1, b2∈{1, 2, …, b}}, 即, 岸桥从贝b1移到贝b2所需时间为|b1-b2|·t, 船贝从左到右依次递增编号

(2) 问题数据

pi:执行任务i所需时间; li:任务i的贝号; rk:岸桥k的就绪时间; l0k:岸桥k的初始贝号; b:最大贝号; t:岸桥行走一个贝所用时间; s:岸桥之间的安全距离; M:一个足够大的正整数。

(3) 决策变量

①tasksi:

任务集合中第i个装卸任务, 定义为区间变量类型, 它具有起始时间、终止时间和大小 (通常情况下, 区间变量的大小等于任务终止时间-任务起始时间) 等内在特征, 在求解结果中将得出其起始时间、终止时间和大小。在OPL中定义形式如下:

dvar interval tasks [i in 1..n] in 0..M size p[i];

该定义说明任务起止时间在0到一个极大的数M之间, 且第i个任务的大小为pi.

②qcTaskskj:

区间变量, 表示岸桥k执行任务j的活动, 其具体定义形式如下:

dvar interval qcTasks[k in 1..q][j in 1..n] optional in r[k]+ftoi (abs (l0[k]-l[j]) ) * t..M;

optional关键词说明岸桥k并不是要执行所有任务的 (即某任务并不一定出现在最终解中) , 具体执行哪些任务将在后面的约束中确定, 该区间变量用于为岸桥分配任务。 qcTaskskj的起始时间最小为:岸桥k的最早可用时间rk+岸桥k从其初始所在贝行走到任务j所在贝的行走时间|l0k-lj|·t, 式中abs函数求绝对值, ftoi函数将浮点数转换成整数。变量范围下界的定义实际上隐含了一个约束情景, 即若岸桥初始位置所对应的船贝并无任务, 它需要消耗时间走到分配给其的第一个任务所在贝处。

③qcsk:

定义为顺序 (sequence) 变量类型, 顺序变量定义在一组区间变量上, 表示这组区间变量的排序。qcsk的定义如下:

dvar sequence qcs[k in 1..q] in all (i in 1..n ) qcTasks[k][i] types all (i in 1..n) l[i];

顺序变量qcsk定义在由岸桥k执行的所有任务上, 其解将是这一组任务的一个排序。还指定了这组任务中, 每个任务都具有一个类型, 用该任务所在的贝号表示。这个任务类型要被后面的约束 (3) 使用。

④qcranei:

整数决策变量, 表示给任务i分配的岸桥编号, 是为方便建模而引入的辅助决策变量。

2.3 约束

(1) 任务顺序关系约束

约束 (1) 确保具有顺序关系的任务在时间上满足顺序要求。

$e n d B e f o r e S t a r t (t a s k s_{i}, t a s k s_{j}), \forall (i, j) \in Ρ h i (1)$

(2) 任务分配和时间约束

约束 (2) 使用了alternative (选择) 约束和操作符all, 其作用是保证如果第i个任务tasksi出现在最终调度中, 那么执行任务i的q个岸桥任务qcTaskski (k∈Q) 中, 只能有一个出现在最终调度中 (即最终执行) , 且该任务的起止时间与tasksi的完全相同 (即, 二者等价) , 这个样就可以保证每个任务都仅被指定给一个岸桥。

$a l t e r n a t i v e (t a s k s_{i}, \underset{k \in Q}{a l l} q c Τ a s k s_{k i}), \forall i \in Τ a s k s (2)$

约束 (3) 利用noOverlap约束保证顺序变量qcsk中的任务在时间上不重叠, 且相继任务间至少要消耗一个过渡时间transitionTimes (即移动时间) , 求解器会将相继任务的类型 (即贝号) 与transitionTimes集合中的元组匹配, 得到相应的过渡时间 (移动时间) 。

$n o Ο v e r l a p (q c s_{k}, t r a n s i t i o n Τ i m e s), \forall k \in Q (3)$

(3) 岸桥冲突约束

岸桥调度中要保证岸桥不能互相跨越, 并且要保持一个安全间距, 对该约束的详细解释参考Bierwirth和Meisel (2009) [1]。

在CP模型中, 有多种方法表达该约束, 效率差别很大, 例如可以用约束 (4a) 表达岸桥冲突约束, 其中的M是足够大的整数。约束 (4a) 的意思是, 假设两个岸桥v<w, 即v在w左边, 分别执行两个任务i和j (任务所在贝满足关系li>lj- (s+1) ·|v-w|) , 则j在i前, 或者i在j前的时间间隔至少为 (li-lj+ (s+1) ·|v-w|) ·t才能保证两个岸桥不互相跨越, 且距离大于最小安全间距, 由于这个约束考虑了j在i前, 或者i在j前两种情况, 故允许岸桥双向移动。endOf函数返回任务的结束时间, 若该任务未出现在最终调度中, 则返回第二个参数值。startOf函数与此类似。约束 (4a) 虽然易于理解, 但因为求解器要迭代v、w、i、j四个量, 执行效率不高。

$\begin{array}{l} e n d Ο f (q c Τ a s k s_{w j}, 0) + (l_{i} - l_{j} + (s + 1) \cdot | v - w |) \cdot t \\ \leq s t a r t Ο f (q c Τ a s k s_{v i}, Μ) | | \\ e n d Ο f (q c Τ a s k s_{v i}, 0) + (l_{i} - l_{j} + (s + 1) \cdot | v - w |) \cdot t \\ \leq s t a r t Ο f (q c Τ a s k s_{w j}, Μ), \\ \forall v, w \in Q | v < w, \\ \forall i, j \in Τ a s k s | i \neq j \land l_{i} > l_{j} - (s + 1) \cdot | v - w | \end{array} (4 a)$

为提高效率, 重新设计一个约束 (4b) 来表达胺桥冲突约束, 经实验验证, 该约束是求解效率更高, 也是本文最终采用的表达岸桥冲突的约束形式。该约束中div是整除操作符, maxl是求最大值的函数。

$\begin{array}{l} q c r a n e_{j} - q c r a n e_{i} > m a x l (0, (l_{j} - l_{i}) d i v (s + 1)) = > \\ e n d Ο f (t a s k s_{j}) + (l_{i} - l_{j} + (s + 1) \cdot (q c r a n e_{j} - q c r a n e_{i})) \cdot t \\ \leq s t a r t Ο f (t a s k s_{i}) | | \\ e n d Ο f (t a s k s_{i}) + (l_{i} - l_{j} + (s + 1) \cdot (q c r a n e_{j} - q c r a n e_{i})) \cdot t \\ \leq s t a r t Ο f (t a s k s_{j}), \\ \forall i, j \in Τ a s k s | i \neq j \end{array} (4 b)$

(4) 不能同时执行的任务约束

约束 (5) 利用函数aqpend将两个任务连成一个顺序变量, 然后利用noOverlap约束保证它们时间上不重叠。

$n o Ο v e r l a p (a p p e n d (t a s k s_{i}, t a s k s_{j}), \forall (i, j) \in Ρ s i (5)$

(5) 设置qcranei

约束 (6) 将执行任务i的岸桥的编号设置给qcranei, 即为任务i分配的岸桥号。如果某任务不出现在最终调度中, 则presenceOf返回0, 否则返回1。

$\begin{array}{l} q c r a n e_{i} = = \min_{k \in Q} (k + (q + 1) (1 - p r e s e n c e Ο f (q c Τ a s k s_{k i}))), \\ \forall i \in Τ a s k s (6) \end{array}$

(6) 累积约束

以上约束已经能够完全定义我们的问题, 但为了提高CP引擎的搜索能力, 这里引入一个累积约束 (7) , 它表示任一时刻所有任务消耗的资源 (这里是岸桥) 总数不得超过q. pulse是OPL中的基本累积函数。当某区间变量qcTaskski表示的活动开始时, pulse函数会导致资源消耗加1, 当该活动结束时, pulse函数会导致资源消耗减1。

$\sum_{k \in Q} \sum_{i \in Τ a s k s} p u l s e (q c Τ a s k s_{k i}, 1) \leq q (7)$

2.4 搜索策略

为了进一步提升ILOG CP优化器的搜索效率, 可利用ILOG CP提供的搜索阶段 (search phase) 机制来干预其搜索策略。通过搜索阶段指定CP优化器固定区间变量的搜索次序, 不同的次序可能会对搜索效率产生显著影响。我们通过实验尝试了几个不同的次序, 发现先固定qcTasks数组, 再固定tasks数组的搜索效率最高。搜索阶段用ILOG脚本定义如下:

cp.setSearchPhases (cp.factory.searchPhase (qcTasks) , cp.factory.searchPhase (tasks) ) ;

2.5 目标函数

目标函数 (8) 是最小化最迟任务的完成时间 (相当于全部任务都完成的时间) 。

$\min \max_{i \in Τ a s k s} {e n d o f (t a s k_{i})} (8)$

3 数值实验

3.1 性能比较

本节比较CP方法与其他方法的求解性能。测试实例由Kim和Park (2004) [4]给出, 包括90个实例, 编号为k13至k102, 这些实例按照规模分成9个实例集 (A到I) , 由于绝大多数文献都测试了A到D实例集, 所以本文也测试A到D实例集。各实例集包含的实例编号、任务数和岸桥数见表1。所有实例的岸桥最早可用时间rk都设为0, 岸桥行走速度t为1个时间单位/贝, 岸桥安全间距s为1贝。

参与比较的算法包括:Kim和Park (2004) [4]的分支定界算法 (B&B) 和贪婪随机适应性搜索算法 (GRASP) 、Moccia等 (2006) [5]的分支割算法 (B&C) 、 Sammarra等 (2007) [6]的禁忌搜索算法 (TS) 、Chung和Choy (2012) [7]的遗传算法 (GA) 、Bierwirth和Meisel (2009) [1]的单向调度启发式算法 (UDS) 、和本文提出的双向调度CP方法。CP模型求解时间限制为1小时。求解结果见表2。表中GA、UDS、CP列数据由笔者整理, 其它列数据由Bierwirth和Meisel (2009) [1]整理记录。

各算法列记录的是该算法找到的最佳值 (fbest) 与实际最优值 (fopt或下界LB) 的相对误差 (RE:relative error) , 即RE = (fbest-fopt) /fopt×100。实例k22和k42在Bierwirth和Meisel (2009) [1]以前的MIP模型中, 没有正确检测出岸桥冲突约束, 因此表中其对应的fopt值稍偏低 (实际分别为540、573) , UDS算法和CP方法在k22、k42得到了实际最优解, 其误差是相对于实际fopu值计算的。

+ 表示下界。 * 此处原始文献中的目标值小于实际最优目标值6.19%, 故原始文献数据应有误。

观察表2, 可以看到, 对最小规模的实例集A, 除了B&B和GRASP外, CP和其它算法都得到了全局最优解。对实例集B, 只有B&C、UDS和CP得到了全局最优解。对实例集C, CP显著优于B&B、GRASP、TS和GA, 略差于B&C和UDS. 对最大的实例集D, CP显著优于B&B、GRASP、B&B、TS和GA, 略差于UDS. 从总体上看, CP方法的结果优于多数算法, 略差于UDS, 但差距非常小。

表3列出了各算法求解每个实例集的平均计算时间 (分钟) 和计算平台, 总体上看CP快于B&B、B&C、GRASP和TS, 慢于GRASP、GA、和UDS, 考虑到CP的求解质量好于GRASP、GA, 且平均求解时间在一刻钟以内, 效果还是不错的。在求解质量几乎相同的情况下, CP慢于UDS的一个原因是UDS是在单向调度的解空间中搜索, 而CP是在双向调度的解空间中搜索, 其搜索空间要大得多。

3.2 双向调度

上述实验结果虽然展示了CP求解质量好于多数算法, 但并未显示出CP模型的可求双向调度解的优势 (因为CP解的质量仍然略差于UDS) , 这可能是因为在测试实例中, 最优解存在于单向调度中, 也可能是因为CP没有搜索到质量更好的双向调度解。目前学界还不清楚具有双向调度最优解的实例特征, 因此无法批量生成具有双向调度最优解的测试实例数据, 从而无法系统比较双向调度下各种算法的质量。但仍然能够构造一些简单的双向调度实例数据, 来检验CP模型是否能够求得双向调度最优解。这里构造一个具有双向调度最优解的实例:有两台岸桥, 岸桥安全间距1个贝, 岸桥移动速度为1, 岸桥就绪时间都为0, 岸桥初始位置在1、3贝。船上有5个贝, 6个装卸任务, 各任务的处理时间分别为11、11、11、22、11、22, 各任务所在贝分别为1、1、3、3、3、5, 具有顺序关系的任务对集合Phi={<1, 2>, <3, 4>, <4, 5>}, 不能同时执行的任务对集合Psi={<1, 2>, <3, 4>, <4, 5>}。

若用UDS算法, 仅能求得单向调度最优解如图1a所示, 目标值为57;而用本文提出的CP方法可以求得双向调度最优解如图1b所示, 目标值为48。图1中横坐标为时间, 纵坐标为贝号, 图中符号如[2][6]表示岸桥2执行任务6的装卸活动。

4 结论

本文构建了一个可求得双向调度解的约束模型用于求解QCSP问题, 数值实验显示其求解质量和效率高于多数现有启发式算法, 略低于UDS算法。但本模型具有求解岸桥双向调度问题的能力, 并在小实例中进行了测试验证。未来需要明确存在双向调度最优解的实例特征, 以便批量生成大型的、具有双向调度最优解的实例对CP方法的性能进行进一步测试。

约束二次规划问题篇5

模拟退火算法来源于固体退火原理, 它是由Kirkpatrick最早提出, 基于Monte-Carlo迭代求解策略的一种随机寻优算法, 其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发, 伴随温度参数的不断下降, 结合概率突跳特性在解空间中随机寻找目标函数的全局最优解, 即在局部最优解能概率性地跳出并最终趋于全局最优。

二次规划是非线性规划中比较简单的一类, 很多方面的实际问题都可以抽象成二次规划的模型去求解, 例如在运筹学中, 它被广泛用于经济调度、合理分配、计划决策等问题。传统解决二次规划问题的方法过程比较复杂, 计算时间较长, 使其在大范围优化中的应用受到限制。文中采用模拟退火算法求解二次规划问题, 具有计算速度快和精度高的特点。

2 设计思路

如图1所示, 考虑如下二次规划问题, 求解下列函数的最小值, xx+yy+xy-4x-6y, 约束条件为x+y=5, 并且x, y都大于零。模拟退火算法的实现过程如下:

(1) 初始化:初始温度Temperature, 初始解状态PreX, 对于每个温度迭代次数MarkovLength。

(2) 对于i=0…MarkovLength, 求解第 (3) 至第 (6) 步。

(3) 产生新解NextX。

(4) 计算增量△=GetFitness (PreX) -Get Fitness (NextX) , 其中GetFitness (doublex) 为评价函数。

(5) 如果△大于零, 则接受NextX作为新的当前解, 否则以概率exp (-△/Temperature) 接受NextX为新的当前解。

(6) 温度逐渐减小, 转到第 (2) 步。

(7) 如果满足终止条件则输出最优解, 结束程序。

网络资料上针对模拟退火算法有很详细的介绍, 但是缺乏翔实的代码实现, 对于初学者很难理解算法和实现, 下文给出基于C#语言的具体实现代码, 读者可以借鉴使用。

3 实现代码

3.1 窗口类代码