调度

管晓宏①, 翟桥柱①*, 冯泳翰②, 高峰①

① 西安交通大学系统工程研究所, 机械制造系统工程国家重点实验室, 智能网络与网络安全教育部重点实验室, 西安 710049; ② 西安交通大学理学院, 西安 710049 * E-mail: qzzhai@sei.xjtu.edu.cn

收稿日期: ; 接受日期:

国家自然科学基金(批准号: 60704033, 60736027)、国家高技术研究发展计划(“863”计划)(批准号: 2007AA04Z1)和教育部新世纪优秀人才计划(批准号: NCET-08-0432)资助项目

摘要 “即时消费”类生产制造系统的优化调度具有重要学术和应用价值. 满足此类系统对产量的实时需求, 考虑调度计划的可实现性具有挑战性. 如何得到精确满足累积产量实时需求的最优调度目前尚无系统方法, 迫切需要研究. 本文建立了含积分约束的生产制造系统优化调度新模型. 通过对生产量变化率约束的深入分析, 证明了该类优化问题等价于光滑非线性规划问题. 生产设备在各时段的产量上下界可表述为时段初、末时刻瞬时生产率的二元函数, 且为精确可达的上下界. 本文结合梯度映射的单调性, 证明了上下界函数的凸性(凹性), 在生产成本为凸函数时, 进一步证明了此类优化调度问题等价于凸规划问题. 本文以上述分析为基础, 针对含积分约束的生产制造系统优化调度问题, 提出了两阶段数值求解方法, 在许多情况下可以迅速获得调度问题的全局最优解. 新模型和相应求解方法克服了生产量变化率约束带来的困难, 获得了精确满足累积产量实时需求的最优调度. 本文同时以电力生产优化调度问题为例, 进行数值求解, 并对结果进行了讨论, 验证了新模型和相应方法的有效性.

关键词 生产调度 积分约束 最优控制 凸规划

生产调度是生产制造系统运行最重要的任务之一. 通常, 调度的主要目的是合理调配资源和设备, 在满足需求、确保安全等运行约束的前提下降低成本、节约资源、减少污染. 优化生产调度能够在不改变基本工艺过程、不更动装备的情况下实现节能降耗, 降低成本, 获得巨大的社会经济效益, 在当前能源、环境问题日趋严峻的形势下, 具有特别重大的意义. 因此, 生产调度理论与方法长期以来一直是一个活跃的研究领域, 受到众多学者关注[1~6].

在生产制造系统中, 有一类生产所谓不可存储产品, 即“即时消费”或“鲜活”(Perishable)类产品的系统[1], 其显著特点是产品必须立即消费而不能存储, 最典型的例子是电力生产. 在电力生产调度问题中[7,8],

每个瞬时生产的电量必须与系统负载需求保持平衡. 其它类似的系统包括管道输送天然气生产和许多服务业的即时服务等. 这类系统的优化调度问题不但会有实时系统需求等, 还可能有十分复杂的离散和连续的动态运行约束, 几乎不可能按连续时间求解. 目前广泛采用的建模方式是将时间离散化, 在一个调度时段内, 假设系统需求为恒定值, 用生产设备的平均生产量满足这一时段恒定的需求. 将平均生产量作为优化决策变量, 求解离散时间点的生产量, 可将调度问题由连续时间最优控制问题转化为数学规划问题, 从而大大降低问题的复杂性[2,7~10].

时间离散化的生产制造系统优化调度问题一般模型为整数规划或混合规划问题, 国内外研究者对

此类问题进行了大量研究[2~19], 取得了许多重要成果. 基于Lagrange松弛的优化方法是解决此类问题最有效的方法之一[3,5,9,11]. 近年来, 随着通用整数规划或混合规划算法效率和计算机性能的大幅度提高, 通用离散和混合优化方法求解此类生产制造系统优化调度问题也取得了很大成功[10,12~14].

由于满足“即时消费”类生产制造系统的实时需求是对调度计划的基本要求, 按离散时间模型获得的调度方案的可实现性就十分重要. 许多生产设备产量的变化率受物理, 如发电机组出力的爬升率有限, 不可能每个时间段起点突变, 完全按上述离散时间模型得到的调度计划不可能操作实现, 也不一定有必要. 实际上很多情况下, 我们只要求这类生产系统在一个时段内的累积产量, 即生产率对时间的积分精确等于实时需求. 例如电力生产系统的生产率是出力(功率), 我们只要求系统在一个时段内交付的能量(功率的积分)满足实时需求. 然而, 在复杂调度问题中“不多不少”精确满足累积产量的实时需求绝非易事. 我们在前期工作中详细阐述了电力生产中离散时间最优调度即使满足出力爬升约束, 也可能存在能量不可交付性, 并给出了判定能量可交付性或调度计划可实现性的充分必要条件[16]. 然而, 如何取得精确满足累积产量实时需求的最优调度目前还没有答案, 迫切需要研究.

本文建立了含积分约束的生产制造系统优化调度新模型. 通过对生产量变化率约束的深入分析, 证明了该类优化问题等价于光滑的非线性规划问题. 生产设备在各时段的产量上下界可表述为时段初、末时刻瞬时生产率的二元函数, 且为精确可达的上下界. 本文结合梯度映射的单调性, 证明了上下界函数的凸性(凹性), 在生产成本为凸函数时进一步证明了此类优化调度问题等价于凸规划问题. 本文以上述分析为基础, 针对含积分约束的生产制造系统优化调度问题, 提出了两阶段数值求解方法, 在许多情况下可以迅速获得全局最优解. 新模型和相应求解方法克服了生产量变化率约束带来的困难, 获得了精确满足累积产量实时需求的最优调度. 本文以电力生产优化调度问题为例, 进行了数值求解, 并对结果进行了讨论, 进而验证了新模型和相应方法的有效性.

应该指出, 本文提出的模型和求解方法可以推广到求解有库存或存储的生产系统调度问题[2]. 限于篇幅, 此类问题未在本文中讨论.

本文结构安排如下: 第1节用一个简单例子说明了目前文献中通用的建模方式存在的缺陷; 第2节给出了按连续时间方式表述生产率及产量积分约束的精确调度模型; 第3节对模型的约束结构特征, 尤其是积分约束的结构特征进行了深入分析, 得到了一些重要理论结果; 第4节证明了新模型等价于一个可行域为闭凸集的光滑非线性规划问题, 特别当生产成本为凸函数时进一步等价于一个简单的凸规划问题, 提出了两阶段数值求解框架; 第5节以电力生产优化调度问题为例验证了本文提出的有关理论和方法; 第6节对全文进行了总结.

1 积分约束的必要性举例

我们以电力生产调度的简单例子, 说明积分约束的必要性.

某台发电机组最小发电功率为150 MW, 最大发电功率为450 MW, 最大爬升速率为6 MW/min, 假定第1 h内的发电功率恒定为g(1) = 150 MW, 现确定第2 h的发电计划.

如果按传统模型, g(1), g(2)表示第1和2 h内机组的平均发电功率(单位: MW). 如每个离散调度时段长度为1 h, g(1), g(2)在数值上等于机组在对应时段的发电量或产出的能量(单位: MWh), 并满足以下约束:

150≤g(2)≤450,g(2)g(1)≤660360.

显然在满足上述约束的情况下, g(2)的最大取值为g(2) = 450 MW. 然而, 容易说明机组在第2 h内的最大发电量或能量远远无法达到450 MWh, 见图1所示.

图1 机组发电功率与电量

因机组在第1 h末发电功率为150 MW, 对应于图

1中A点, 由于爬升速率有限, 发电功率不可能在第2 h初突变, 即便按最大爬升速率上升, 至多在1 h 50 min处达到最大发电功率450 MW, 对应图中C点, 此后保持不变. 因此机组在第2 h内最大可实现的发电量为五边形 ACDFE 的面积, 对应为 325 MWh, 远低

于450 MWh(对应于矩形 BDFE 的面积). 因此, 如果要进行电量的调度, 应该考虑发电功率的积分和相应的约束.

2 含积分约束生产制造系统优化调度模型

考虑一个生产制造系统有I台生产设备, 生产过

程需消耗J种原料、调度周期为K时段. 本文采用的变量和符号定义及说明如下:

i 设备编号, i1,2,,I;

j 生产原料编号, j1,2,,J; k

调度时段编号,

k1,2,,K; 

每个调度时段的时间长度;

aj,i

第i台设备生产每单位的产品, 需消耗的第j种原料量;

Rj,k 第k时段允许消耗的第j种原料总量; t 连续时间变量, 从0时刻开始计时, 调度周期总时间长度为K;

gi(t) 设备i在t时刻的生产率, 是时间的连续函数;

ui(t) 设备i在t时刻的生产率变化率, ui(t)是分段连续函数;

pi(k) 设备i在第k时段的生产量; i

设备i的生产率变化率上限; gi,gi

设备i的最大、最小生产率; Ci(pi(k)) 设备i在第k时段生产量为pi(k)时的生

产成本, 假定成本仅与生产量有关;

D(k)

第k时段的总需求或合同供货量.

本文考虑的生产调度可描述为以下优化问题. 目标: 生产成本最小 IK

uminJi(t),gi(t),pi(k)Ci(pi(k)).

(1)

i1k1约束条件:

(a) 即时产需平衡:

I

pi(k)D(k), k1,2,,K;

(2)

i1(b) 原料: I

aj,ipi(k)≤Rj,k,

(3)

i1其中j1,2,,J, k1,2,,K. 约束(3)具有通用性, 可以表示许多形式的约束, 除原料约束外, 还可表示电力生产调度中的系统备用约束、传输容量约束(均可视为广义的原料约束)等.

(c) 产量-生产率关系:

pi(k)k(k1)gi(t)dt,

(4)

其中i1,2,,I, k1,2,,K.

(d) 生产率-生产率变化率关系:

ggti(t)i(0)0ui()d,t[0,K],

(5)

其中i1,2,,I. 在本文考虑的调度问题中, ui(t)为决策变量, gi(t)和pi(k)为状态变量.

(e) 生产率变化率界约束(生产率爬升约束):

i≤ui(t)≤i,t0,K,

(6)

其中i1,2,,I. (f) 生产率界约束:

gi≤gi(t)≤gi,t0,K, (7)

其中i1,2,,I.

(g) 调度初始及终止时刻生产率约束:

gi(0)gi,0,gi(K)gi,K,

(8)

其中gi,0和gi,K是给定常数. 该组约束中的后一个可选, 某些系统可能不需要.

由(1)~(8)式描述的具有实际背景的生产制造系统优化调度问题十分复杂, 既包括了离散时间约束, 也包括了连续时间积分关系, 难以直接应用现有的非线性规划或最优控制理论求解, 本文通过深入研究约束(4)~(7)及目标函数(1)的结构特征, 找到了将此类调度问题转化为光滑凸规划问题的系统方法.

受篇幅, 模型(1)~(8)中未引入描述设备开关机状态的离散决策变量和约束. 作者已结合本文的理论结果和前期研究成果[9], 在上述模型中考虑了离散决策变量和约束, 将在以后报告.

3 问题转化为光滑凸规划的理论基础

本节给出将问题(1)~(8)转化为光滑凸规划的理论基础. 为使表述简洁, 引入下述符号:

gi,kgi(k),k0,1,2,,K,i1,2,,I,

(9)

表示设备i在k时刻的瞬时生产率.

首先, 可以证明单台生产设备在各时段的产量上下限可表示为该设备在相应时段初、末时刻生产率的二元函数.

定理 1. 设 pi(k), gi(t), ui(t)满足约束(4)~(7), 则有下述结论成立:

(i)

gi,kgi,k1≤i;

(10)

(ii) Pi(gi,k1,gi,k)≤pi(k)≤Pi(gi,k1,gi,k). (11) (11)式中的产量上下限二元函数Pi(,)和Pi(,)解析表达式如下:

Pi(gi,k1,gi,k)(gi,k1g)2(gi,kg)2iigi,

2iifgi,k1gi,k2gii, (gi,k1gi,k)2(gg4i,k1i,k)i2,i24ifgi,k1gi,k≥2gii,(12)Pi(gi,k1,gi,k)(gigi,k1)2(gigi,k)22gi,

iifgi,k1gi,k≥2gii, (gi,k1gi,k)2(gi,k1g)i2,4i,ki24ifgi,k1gi,k2gii.(13)(iii) (11)式给出的pi(k)的上下界是可以达到的, 即分别存在适当的生产率函数gi(t)和生产率变化率函数ui(t)使得(11)式中的两个不等式成为等式.

证明:

(i) 在(5)式中分别令t(k1)和tk可得:

gi,k1gi((k1))g(k1)i(0) 0ui()d,(14)

gi,kgi(k)gk

i(0)0ui()d,(14)式中的两式相减可得:

gi,kgi,k1k(k1)ui()d,

(15)

再根据(6)式, 结合(15)式即可得(10)式, 结论(i)因此成立.

(ii) 当t[(k1),k]时, 对比约束(5)和(6)式, 可得如下结论:

gi(t)gi,k1t(k1)ui()d

≥gti,k1(k1)id

gi,k1i(t(k1)),(16)

gi(t)gi,k()d

ktui

≥gi,kktidsgki(kt),(17)

因此结合(7)式可知, 当t[(k1),k]时下式成立:

gi(t)≥gimin,k(t)

maxgi,k1i(t(k1)),gi,ki(kt),gi, (18) 上式右端给出的gimin,k(t)实际上是在给定第k时段初、末时刻生产率为gi,k1和gi,k的前提下, 在第k时段设备i的最低可达生产率曲线. 根据产量-生产率积分关系约束(4)式, 从而得到:

pktki(k)(k1)gi(t)d(k1)gimin,k(t)dt

(gi,k1gi)2(gi,kgi)2gi2,i

ifgi,k1gi,k2gii, (gi,k1gi,k)2(gg)i24i2i,k1i,k4,ifgi,k1gi,k2gii.(19)(19)式的结果基于(18)式, 通过图2中最左边一列的两幅图(最小可达生产率曲线)能够更好说明. 图2中的Case 1和Case 2分别对应于(19)式中的第一种情况和第二种情况.

基于完全类似的推理分析, 可以得到当t [(k1),k]时有下式成立:

图2 设备i在时段k的最低/最高可达生产率曲线及其对应的生产爬升率

gi(t)≤gimax,k(t)

mingi,k1i(t(k1)), gi,ki(kt),gi, (20) 上式右端给出的gimax,k(t)实际上是在给定第k时段初、末时刻设备生产率为gi,k1和gi,k的前提下, 在第k时段设备i的最高可达生产率曲线. 根据产量-生产率积分关系式(4), 得到:

将变为等号. 本文仅指出对应于图2中与gimin,k(t)对应

的Case 1情况下, ui(t)与gi(t)的对应关系, 其余情况完全类似.

在所考虑的情况下有下式成立:

gi,k1gi,k2gii,

(22)

此时有: t[(k1),k]时

pi(k)k(k1)gi(t)dt≤k(k1)gimax,k(t)dt

(gigi,k1)2(gigi,k)2gi,2iifgi,k1gi,k≥2gii,

2i2(gi,k1gi,k)(gg),i,k1i,k424iifgi,k1gi,k2gii,(21)gi,k1i(t(k1)),if(k1)≤t≤tk,1,gimin(t)gi,iftk,1≤t≤tk,2,,kgi,ki(kt),iftk,2≤t≤k,(23)

其中的tk,1, tk,2由下式得到:

tk,1(k1)(gi,k1gi)/i, tk(gg)/,i,kiik,2(24)

不难验证(22)式满足时, (24)中给出的两个时刻tk,1, tk,2均位于第k时段. 与gimin,k(t)对应的ui(t)为: t[(k1), k]时

(21)式的积分基于(20)式, 通过图2中的中间一列两

幅图(最高可达生产率曲线)能够更好说明. 图2中的Case 1和Case 2分别对应于(21)式中的第一种情况和第二种情况. 结论(ii)因此成立.

(iii) 显然, 在(ii)的证明中, 如果当t[(k1),k]时取gi(t)gimin,k(t), 则(19)式中的不等号将变为等号; 同理, 如果取gi(t)gimax,k(t), 则(21)式中的不等号

i,if(k1)≤ttk,1, ui(t)0,iftk,1≤ttk,2,i,iftk,2≤t≤k.(25)

(23)~(25)式的结果可以用图2中第3列的两幅图

来解释. 注意(25)式中给出的生产率变化率(控制量)曲线呈阶梯状. 定理1至此全部证完.

定理1是将具有积分约束的调度问题(1)~(8)转化为非线性规划问题的基础. 定理2进一步指出当(11)~(13)式中给出已知时段初、末时刻生产率时, 时段内产量上下界函数为凸或凹.

定理2. 定理1中给出的二元函数Pi(,)是光

滑的(即连续可导)凸函数, Pi(,)是光滑的凹函数.

为证明定理2, 需要两个引理.

引理1. 假设f1:RnR, f2:RnR均为连

续可导函数, Hx|Tx,xRn是Rn中的一个超平面, 其中Rn, R为给定向量和实数. 如果f1(x)f2(x)且f1(x)f2(x) 对所有xH成立, 即两个函数在超平面H上光滑衔接, 则下述函数在Rn上连续可导.

f(x)f1(x),ifTx≤, f2(x),ifTx.(26)

证明: 根据导数和光滑函数的定义, 引理1的

结论是显然的.

引理2. 所有条件同引理1, 再假定f1, f2在Rn

中均为凸函数, 则(26)中给出的函数f(x)也是Rn中的凸函数.

证明: 首先, 根据已知条件可得下式成立:

f1(x),ifT

f(x)x≤,f2(x),ifTx≥,T (27)

f(x)f1(x),ifx≤,f2(x),ifTx≥,根据光滑函数成为凸函数的充分必要条件, f1和f2均为Rn中的单调映射[20], 即对任意x1, x2Rn有下式成立:

T

f1(x1)f1(x2)(x1x2)0,

(28)

f2(x1)f2(x2)T(x1x2)0,为证f是凸函数, 只需证明f是Rn中的单调映射, 即对任意x1,x2Rn,

f(x1)f(x2)T(x1x2)≥0,

(29)

(29)式的证明可基于对x1, x2和超平面H间的所有可能位置关系进行分析来完成. 所有可能位置关系全部在表1中列出.

表1 x1, x2和超平面H间的所有可能位置关系

关系分类

Tx2 Tx2 Tx2 Tx1 Case 1 Case 1 Case 2 Tx1

Case 1 Case 1 Case 1 Tx1

Case 2

Case 1

Case 1

在表1中, 所有可能的位置关系被分为两大类, 第一类为:

Case 1. Tx1≤,Tx2≤或Tx1≥,Tx2≥, (30) 在这种类型中, 两个点x1, x2位于超平面H的同一侧, 因此由(27)式可知

f(x1)f1(x1),f(x2)f1(x2),

ifTx1≤,Tx2≤,f(x1)f(x1),f(x2)f(x2), 22ifTx1≥,Tx2≥,(31)

再结合(28)式即可知(29)式成立.

第二类为: Case 2.

Tx1,Tx2或Tx1,Tx2, (32)

在这种类型中, 两个点x1, x2位于超平面H的不同侧.

基于问题对称性并不失一般性, 仅考虑Tx1>, Tx2>发生的情况. 此时由(27)式可知:

f(x1)f1(x1),f(x2)f2(x2),

(33) Tx1令 Tx2Tx1,x0x1(x2x1), (34)

则有

01,Tx0.

(35)

再根据f1和f2的凸性可知:

T

f1(x1)f1(x0)(x1x0)≥0, (36)

f2(x0)f022(x2)T(xx)≥0.由(34)式可得:

x1x0(x1x2),x0x2(1)(x1x2),

将上式代入(36)式并注意到01则可以得到:

T

f1(x1)f1(x0)(x1x2)≥0,f02T (37)

122(x)f2(x)(xx)≥0,因为x0H(见(35)式)且根据引理的条件f1和f2在H上一致, 从而有f1(x0)f2(x0). 因此将(37)式中两式相加即得:

f12T121(x)f2(x)(xx)≥0,

(38)

此即(29)式. 引理2至此证完.

以上述两个引理为基础, 下面给出定理2的证明. 证明: 我们仅证明Pi(,)是连续可导的凸函数, 用完全类似的方法可以证明Pi(,)是连续可导的凹函数(证明Pi(,)为连续可导凸函数即可). 令: f1(gi,k1gi,k1,gi,k)(gi)2(gi,kg2i)2gii, (39)

(gi,k1gi,k)2f(g2i,k1,gi,k)4(gi,k1gi,k)i2,i24(40) 显然f1, f2均为R2上的连续可导下凸函数, 现在考虑R2中的下述超平面(直线):

H(gi,k1,gi,k)|gi,k1gi,k2gii, (41) 根据(39)和(40)式可得:

gi,k1gif1(gi,k1,gi,k)igi,kg,i igi,kg (42)

1i,kif2i2(gi,k1,gi,k)gi,kgi,,k1i2i由此可得:

ifgi,k1gi,k2gii,thenf1(gi,k1,gi,k)f2(gi,k1,gi,k), (43)

即f1,f2在超平面H上保持一致. 再令:

1(y)f1(y,2giiy),(44)

2(y)f2(y,2g iiy),因此, 根据链导法则及(43)和(44)式可以验证下式成立:

1(y)11f1(y,2giiy)11f2(y,2giiy)(y),

2(45)

其中“”表示向量内积运算. 此外, 容易验证:

1(gii/2)2(gii/2),

(46)

(45)和(46)式结合起来表明两个一元函数1(y)和

2(y)不仅导数始终相等, 而且在某一点上函数值相

同, 因此这两个函数实际上恒等, 即:

1(y)2(y),

从而由(44)式可得:

ifgi,k1gi,k2gii,thenf1(gi,k1,gi,k)f2(gi,k1,gi,k), (47)

上式表明f1, f2在H上保持一致, 再结合(43)式可知f1, f2在H上光滑衔接. 又因为:

f1(gi,k1,gi,k),Pifgi,k1gi,k2gii,i(gi,k1,gi,k) f2(gi,k1,gi,k),ifgi,k1gi,k≥2gii,(48)所以根据引理1和引理2, Pi(,)是连续可导的下凸

函数. 定理2至此全部证完.

4 求解积分约束优化调度问题的系统算法

基于第3节的理论结果, 本节将证明最优控制形式的调度问题(1)~(8)可转化为一个普通的光滑非线性规划问题, 对某些重要应用可进一步转化为凸规划问题.

考虑下列非线性规划问题:

IK gminJi,k,pi(k)Ci(pi(k)),

(49) i1k1 s.t.Ipi(k)D(k);

k1,2,,K,

(50) i1Iaj,ipi(k)≤Rj,k;j1,2,,J, k1,2,,K,

(51) i1

gi,kgi,k1≤i,

(52) Pi(gi,k1,gi,k)≤pi(k)≤Pi(gi,k1,gi,k), (53) gi≤gi,k≤gi,

()

gi,0gi,0,gi,Kgi,K,

(55)

约束(52)~()中下标的变化范围为: i1,2,,I, k1,2,,K. (53)式中产量上下限二元函数Pi(,)和Pi(,)由(12)和(13)式确定; (55)式中的gi,0和

g i,K是给定常数, 对某些系统该式不需要考虑,与(8)

式对应. 我们先将问题(49)~(55)看作的非线性规划问题, 再建立问题(1)~(8)与(49)~(55)的关系.

定理3. 下述两个结论成立.

(ⅰ) 设pi(k), gi(t), ui(t)是问题(1)~(8)的一个解, 则pi(k), gi,k是问题(49)~(55)的一个解, 且对应的目标函数值相同, 其中gi,k按(9)式得到.

(ⅱ) 设pi(k), gi,k是问题(49)~(55)的一个解, 则由此可以得到问题(1)~(8)的一个解pi(k), gi(t), ui(t), 且对应的目标函数值相同. 其中gi(t)与gi,k满足(9)式.

定理3表明求解具有积分约束的优化问题(1)~(8), 可以转化为非线性规划问题(49)~(55). 问题(49)~(55)没有与决策变量ui(t)相直接对应的变量, 除去约束(53)外, 其余约束均为线性约束, 因此问题结构比较简单. 在对结论(ⅱ)证明时, 可由(49)~(55)的一个解构造(1)~(8)的一个解, 且便于实现. 定理3的证明如下.

证明: 对于结论(ⅰ), 根据问题(1)~(8)及问题(49)~(55)的表述形式, 再结合定理1的前两个结论可知约束(52)和(53)成立, 其他约束成立是显然的, 两个问题的目标函数完全相同, 因此结论(ⅰ)成立. 下面证明结论(ⅱ).

设pi(k), gi,k是问题(49)~(55)的一个解, 则由定理1的证明过程(参考(18)和(20)式)及(53)式可知:

gi≤gimin,k(t)≤gimax,k(t)≤gi,P)ki(gi,k1,gi,kmin(k1)gi,k(t)dt≤pi(k)(56) ≤k(k1)gimax,k(t)dtPi(gi,k1,gi,k),另外, 在定理1结论(iii)的证明中, 以gi(t)gimin,k(t)为例指出了如何构造ui(t)的方法. 利用同样的构造方

法, 设t[(k1),k]时, 与gimin,k(t)和gimax,k(t)对应的ui(t)分别为uimin,k(t)和uimax,k(t), 则有下述结论成立:

t[(k1),k]时,

i≤uimin,k(t)≤i,i≤uimax,k(t)≤i, gtimin,k(t)gi,k1(k1)uimin,k()d, (57)

gtimax,k(t)gi,k1(k1)uimax,k()d,令

0, ifPi(gi,k1,gi,k)Pi(gi,k1,gi,k)pi(k)Pi(gi,k1,gi,k)kPi(gi,k1,gi,k)Pi(gi,k1,g, i,k)ifPi(gi,k1,gi,k)Pi(gi,k1,gi,k),(58)由(53)及(58)式可知

0≤k≤1, (59)

接下来, 对k1,2,,K, 按如下方式依次在区间t[(k1),k]上构造ui(t)和gi(t):

ui(t)(1k)uimin,k(t)kuimax,k(t),

(60)

gi(t)gi,k1t(k1)ui()d,

(61)

根据(57)和(60)式, (61)式可改写为: 在区间t

[(k1),k]上

gi(t)(1k)gimin,k(t)kgimax,k(t).

(62)

按以上方式得到t[(k1),k](对应于k=1, 2 ,K的各个区间)上ui(t)和gi(t)的表达式后, 就完

成了ui(t)和gi(t)在[0,K]整个区间上的构造.

对于生产率gi(t)而言, 根据(18)式和(20)式可知:

gimin,k(k)gimax,k(k)gi,k,

(63)

将上式代入(62)式可知, 按gi(t)在区间t[(k1), k]上的定义, 有:

gi(k)kgimin,k(k)(1k)gimax,k(k)gi,k. ()

可以看出生产率虽然按区间逐段定义, 但在相邻区间的端点处可以连续拼接, 上述结果同时表明gi(t)与gi,k满足(9)式.

必须注意, 对ui(t)进行拼接后, 在相邻区间的端点处ui(t)可能不连续, 但对生产率的变化率是允许的, 相当于最优控制中的“梆-梆”控制, 对最终结果并无影响. 顺便指出, 定理1证明中构造的ui(t)本身就是分段连续函数, 存在不连续点.

以上仅解释了如何基于pi(k)和gi,k构造gi(t) 和ui(t), 下面证明如此构造出的gi(t)和ui(t)对应于(1)~(8)的一组可行解, 且目标函数值相同.

由于问题(49)~(55)的目标函数与(1)~(8)的目标函数完全相同, 我们仅需证明pi(k), gi(t), ui(t)满足系统(1)~(8)中的所有约束.

约束(2), (3), (8)直接对应于约束(50), (51), (55), 因此仅需证明约束(4)~(7)成立. 对于约束(4), 由(53)式和k的定义((58)式)可知: pi(k)(1k)Pi(gi,k1,gi,k)kPi(gi,k1,gi,k). (65)

将(56)式中的第二式代入(65)式, 再结合(62)式即可得(4)式.

约束(6)式可由(57)式第一行、(59)和(60)式结合得到. 约束(7)式可由(18), (20), (59), (62)式结合得到. 约束(5)式可基于(61)式对k用归纳法得到. 至此, 定理3全部证完.

定理3表明可以通过求解一个光滑的非线性规划问题来获得原生产调度问题的最优解. 通过对问题(1)~(8)的进一步分析, 我们发现如果原生产调度问题的目标函数为凸函数, 则等价的非线性规划问题(49)~(55)是光滑凸规划问题. 这是非常重要的理论结果, 因为很多实际问题的目标函数为凸函数, 其任意局部最优解也是全局最优解, 能够非常有效求解. 一些经典算法如内点法、序列二次规划、广义简约梯度法等方法都能非常迅速地收敛到全局最优解. 定理4解释了问题(49)~(55)与凸规划的联系.

定理4. 下述两个结论成立.

(ⅰ) 问题(49)~(55)的可行域是一个闭凸集. (ⅱ) 若Ci()(1≤i≤I)均为(光滑)凸函数, 则

(49)~(55)式是一个(光滑)凸规划问题.

证明: 根据问题(49)―(55)的结构特点, 可以发现约束(50)~(55)中, 除约束(53)外, 其余约束均为优化变量的线性等式/不等式约束, 所以如果能证明约束(53)的凸性, 结论(i)即可得证.

将约束(53)改写为:

Pi(gi,k1,gi,k)pi(k)≤0,

(66)

pi(k)Pi(gi,k1,gi,k)≤0,

(67)

根据定理2, Pi(,)是光滑的凸函数, Pi(,)是光滑的凹函数, 从而Pi(,)是光滑的下凸函数, 由此观察约束(66)和(67)可知(66)和(67)的左端均为决策

变量的凸函数, 又因为这两个约束都是“≤0”型约束, 因此问题(49)~(55)的可行域是一个闭凸集.

在结论(i)的基础上, 结论(ii)是显然的. 定理4证毕. 在电力系统生产调度等很多实际问题中, 生产制造成本或费用是凸函数, 如火电机组的煤耗曲线常取为二次凸函数, 即Ci()为凸函数, 此时求解非

线性规划问题(49)~(55)可以得到全局最优解. 按照定理3结论(ii)证明中给出的构造方法, 可以重构出ui(t)和gi(t), 进而得到了具有积分约束的优化问题

(1)~(8)的全局最优解.

基于本节的理论结果, 本文设计了两阶段求解方法, 对问题(1)~(8)进行数值求解. 先采用常规算法求解非线性规划问题(49)~(55), 再依据其解按照定理3结论(ⅱ)证明中给出的方法重构出问题(1)~(8)的解. 由定理4结论(ⅰ), 无论目标函数性质如何, 问题(49)~ (55)的可行域始终为闭凸集, 因此在求解非线性规划时算法将迅速收敛于局部最优解; 特别地, 当生产成本函数为凸函数时, 将收敛于全局最优解. 图3描述了这种两阶段框架.

在上述算法框架下, 第一阶段求解的是一个一般的非线性规划问题(可能为凸规划), 因此任何适合于一般非线性规划的算法都可以用于该阶段问题的求解. 在本文的数值测试中采用的是序列二次规划法[21], 其详细步骤在一般的非线性规划专著中均可找到. 因此这里仅给出第二阶段的详细步骤, 我们将其总结为下述算法.

第二阶段算法.

Step 0. (初始化)记录第一阶段得到的解, 即: pi(k), gi,k(i1,2,,I, k1,2,,K). 然后, 置k=0.

Step 1. 若 k = K, 停止; 否则置 k =k+1继续. Step 2. 根据(12)~(13)式计算Pi(gi,k1,gi,k)和

Pi(gi,k1,gi,k).

Step 3. 按定理3证明中指出的方法在区间

t[(k1),k上计算]uimin,k(t)和uimax,k(t)(阶梯函数).

Step 4. 根据(58)式计算k.

Step 5. 根据(60)―(61)式在区间t[(k1),k]上计算ui(t)和gi(t), 返回Step 1.

图3 两阶段法框架

以上两阶段算法最主要的特点体现在两个方面: 第一, 问题模型准确, 考虑了生产率积分约束, 保证了产量优化调度结果的可实现性, 克服了传统模型的缺陷; 第二, 当生产成本为凸函数时, 得到的是全局最优解.

成本(发电煤耗成本)函数表达式为:

Ci(pi(k))ai(pi(k))2bipi(k)ci,

(68)

其中ai,bi,ci均为非负常数, 因此生产成本为产量的光滑凸函数. 各机组的主要参数及系统在各时段的负荷需求分别在表2和表3中给出.

我们在PⅣ2. 0GHz Windows工作站上应用Matlab 6.5优化工具箱中的序列二次规划函数对问题求解, 时间约为418 s, 获得机组1和机组3的最优生产率曲线如图4所示.

在图4中, 每个圆圈标示出了在相邻时段交界点处的机组瞬时出力. 对比系统负荷需求数据(表3)和机组参数(表2)可以发现, 机组1的生产成本最低, 因此在负荷需求较大时处于满负荷运转; 机组3的生产成本较高, 因此仅在负荷需求较大时发电功率较大, 在负荷需求很低时发电功率处于较低水平. 对于机

5 应用示例

我们将本文的理论结果应用于电力系统生产优化调度(亦称经济分配), 经大量实例测试, 验证了本文提出的求解方法非常有效, 本节简要介绍8机组优化调度案例的测试结果.

考虑一个8台机组的电力系统生产优化调度问题, 调度周期为24 h, 每个时段1 h. 对应于模型(1)~ (8)有: I8,K24h,1h. 标准的离散时间经济分配模型可参见文献[22,23]. 各机组(生产设备)的生产

表2 各机组物理参数

Unit No. (i) 1 2 3 4 5 6 7 8

gi(MW) gi(MW) i (MW/h) 600.6 600.6 247 243 243 240 144 102.3

ai ($/(MWh)2) 0.00031 0.00031 0.00395 0.00395 0.00399 0.00398 0.00712 0.00222

bi ($/MWh) 17.62 17.62 19.5 19.7 19.8 19.7 22.26 27.27

ci ($) 970 970 456 450 445 450 370 665

gi(0) (MW) 300 300 60 55 55 60 30 10

455 455 180 170 162 162 80 65 150 150 25 25 25 25 25 10

表3 系统在各时段的负荷需求

k

D(k) (MWh)

k

D(k) (MWh)

1 765 13 1470

2 805 14 1420

3 875 15 1360

4 915 16 1360

5 976 17 1390

6 1055 18 1450

7 1265 19 1370

8 1270 20 1320

9 1465 21 1350

10 1465 22 1300

11 1683 23 840

12 1688 24 760

图4 机组1(a)和机组3(b)的最优出力曲线

组3, 初始时刻发电功率较高是因为受到初值约束(约束(8)).

在此必须说明, 对于第一阶段得到的pi(k)和gi,k, 可能有无穷多gi(t)和ui(t)使得约束(4)~(9)成立, 即第二阶段可能有无穷多组解. 这种现象已在定理3的证明中有所暗示. 定理3的证明是一种构造性证明, 仅指出了解的存在性和一种最自然的构造方法, 并未讨论解的唯一性. 上节的算法中给出的也是最易于理解和编程实现的一种构造算法. 在实际应用中, 可以根据可能的二级优化目标或其他要求从多组gi(t)和ui(t)中选出最合适的一组作为最终解. 在本节的算例中, 我们基于一种系统化算法从众多的gi(t)和ui(t)中选出了一组最“光滑”的gi(t)及其对应的ui(t)作为最终解, 其物理意义为使机组出力尽量平稳, 在满足系统需求和总成本最低的前提下机组出力爬升尽量小. 由于详细的实现过程与本文主题关系不大, 对此问题我们已另文讨论.

注意到图4中的横轴时间单位为小时, 最优的机组出力曲线实际上相当平稳, 完全满足爬升约束等对机组出力曲线的要求. 算例测试结果表明, 本文提出的方法是有效的, 得到了最优的生产调度计划.

6 结论

在“即时消费”型产品的生产系统优化调度中, 积分约束通常必须考虑. 目前广泛采用的离散时间模型可能导致产量计划不可实现. 本文分析了现有模型的缺陷, 建立了具有积分约束的优化问题模型, 并证明此类问题可转化为光滑的非线性规划问题, 当费用函数为凸(或效益函数为凹)时可进一步等价为光滑凸规划问题. 本文给出了构造原问题等价最优解的系统方法, 并提出了两阶段求解框架. 新模型及其求解算法不仅克服了传统模型下的调度结果不可实现问题, 而且当生产成本为凸函数时可以获得全局最优解. 基于电力系统调度的实例测试表明本文提出的相关理论和算法非常有效.

参考文献

1 Talluri T T, Van Ryzin G J. The Theory and Practice of Revenue Management. Heidelberg: Springer, 2004

2 Bannister C H, Kaye R J. A rapid method for optimization of linear systems with storage. Oper Res, 1991, 39: 220—232 3 Chen H, Chu C, Proth J M. An improvement of the Lagrangian relaxation approach for Job Shop Scheduling: A dynamic

programming method. IEEE Trans Rob Autom, 1998, 14: 786—795

4 Cohen A I, Sherkat V. Optimization-based methods for operations scheduling. Proc IEEE, 1987, 75: 1574—1591 5 王朝晖, 陈皓勋, 胡保生. 用Lagrangian松弛法解化工批处理调度问题. 自动化学报, 1998, 24: 1—8

6 Muiser R F H, Evans L B. An approximated method for the production scheduling of industrial batch process with

parallel units. Comp Chem Eng, 19, 13: 229—238

7 Salam M S, Nor K M, Hamdan A R. Hydrothermal scheduling based Lagrangian relaxation approach to hydrothermal

coordination. IEEE Trans Power Syst, 1998, 13: 226—235

8 Bard J F. Short-term scheduling of thermal-eglectric enerators using Lagrangian relaxation. Oper Res, 1988, 36: 756—

766

9 翟桥柱, 管晓宏, 郭燕, 等. 具有混合动态约束的生产系统优化调度新算法. 自动化学报, 2004, 30: 539—6

10 Sand G, Engell S. Modeling and solving real-time scheduling problems by stochastic integer programming. Comp Chem Eng,

2004, 28: 1087—1103

11 Guan X, Guo S, Zhai Q. The conditions for obtaining feasible solutions to security-constrained unit commitment

problems. IEEE Trans Power Syst, 2005, 20: 1746—1756

12 Fu Y, Shahidehpour M. Fast SCUC for large-scale power systems. IEEE Trans Power Syst, 2007, 22: 2144—2151

13 Silva B, Stursberg O, Krogh B, et al. An assessment of the current status of algorithmic approaches to the verification of

hybrid systems. Proceedings of IEEE Conference on Decision and Control. Orlando: IEEE, 2001, 12: 2867—2874 14 Till J, Engell S, Panek S, et al. Applied hybrid system optimization: An empirical investigation of complexity. Control

Eng Practice, 2004, 12: 1291—1303

15 Ferreira L A F M, Anderson T, Imparato C F, et al. Short-term resource scheduling in multi-area hydrothermal power

systems. Electric Power & Energy Systems, 19, 11: 200—212

16 Guan X, Gao F, Svoboda A J. Energy delivery capacity and generation scheduling in the deregulated electric power

Market. IEEE Trans Power Syst, 2000, 15: 1275—1280

17 Wang C, Shahidehpour S M. Optimal generation scheduling with ramping costs. IEEE Trans Power Syst, 1995, 10: 60—

67

18 余贻鑫, 王东涛. 输电系统动态安全风险评估与优化. 中国科学E辑: 技术科学, 2009, 39: 286—292

19 Si B F, Long J C, Gao Z Y. Optimization model and algorithm for mixed traffic of urban road network with flow

interference. Sci China Ser E-Tech Sci, 2008, 51: 2223—2232

20 Hiriart-Urruty J, Lemarechal C. Fundamentals of Convex Analysis. Heidelberg: Springer, 2001

21 Bazaraa M S, Sherali H D, Shetty C M. Nonlinear Programming: Theory and Algorithms. 2nd ed. New York: John Wiely,

1993

22 Travers D, Kaye R J. Dynamic dispatch by constructive dynamic programming. IEEE Trans Power Syst, 1998, 13: 72—

78

23 Han X S, Gooi H B, Kirschen D S. Dynamic economic dispatch: Feasible and optimal solutions. IEEE Trans Power Syst,

2001, 16: 22—28