A Hybrid MILP and IPM for Dynamic Economic Dispatch with Valve Point Effect

TL;DR

This paper introduces a hybrid MILP and IPM approach to efficiently solve the complex, non-convex dynamic economic dispatch problem with valve-point effects, including transmission loss, improving solution quality.

Contribution

A novel hybrid MILP-IPM method is proposed, combining MILP for initial solution and IPM for refinement, effectively handling non-differentiability and non-convexity in DED-VPE.

Findings

The method effectively solves DED-VPE with valve-point effects.

Simulation results validate the approach's accuracy and efficiency.

The hybrid approach outperforms traditional methods in solution quality.

Abstract

Dynamic economic dispatch with valve-point effect (DED-VPE) is a non-convex and non-differentiable optimization problem which is difficult to solve efficiently. In this paper, a hybrid mixed integer linear programming (MILP) and interior point method (IPM), denoted by MILP-IPM, is proposed to solve such a DED-VPE problem, where the complicated transmission loss is also included. Due to the non-differentiable characteristic of DED-VPE, the classical derivative-based optimization methods can not be used any more. With the help of model reformulation, a differentiable non-linear programming (NLP) formulation which can be directly solved by IPM is derived. However, if the DED-VPE is solved by IPM in a single step, the optimization will easily trap in a poor local optima due to its non-convex and multiple local minima characteristics. To exploit a better solution, an MILP method is required to solve the DED-VPE without transmission loss, yielding a good initial point for IPM to improve the quality of the solution. Simulation results demonstrate the validity and effectiveness of the proposed MILP-IPM in solving DED-VPE.

Tables

Table 1. TABLE I: Summary results for the 5-unit system without loss

Method	Total generation cost ($)			S-time(min)
	Minimum	Average	Maximum
MILP	42563			0.82
MILP-IPM	42524			0.86

Table 2. TABLE II: Outputs (MW) for the 5-unit system without loss

$t$	unit $1$	unit $2$	unit $3$	unit $4$	unit $5$
1	16.7925	98.5398	30.0000	124.9079	139.7598
2	10.0000	98.5398	61.7925	124.9079	139.7598
3	10.0000	98.5398	101.7925	124.9079	139.7598
4	10.0000	98.5398	112.6735	124.9079	183.8788
5	10.0000	80.8990	112.6735	124.9079	229.5196
6	40.0000	100.5398	112.6735	125.2671	229.5196
7	10.0000	98.5398	112.6735	175.2671	229.5196
8	10.0000	91.9911	112.6735	209.8158	229.5196
9	39.4513	98.5398	112.6735	209.8158	229.5196
10	53.4513	98.5398	112.6735	209.8158	229.5196
11	69.4513	98.5398	112.6735	209.8158	229.5196
12	75.0000	112.9911	112.6735	209.8158	229.5196
13	53.4513	98.5398	112.6735	209.8158	229.5196
14	39.4513	98.5398	112.6735	209.8158	229.5196
15	34.0000	98.5398	112.6735	179.2671	229.5196
16	10.0000	98.5398	112.6735	129.2671	229.5196
17	10.0000	80.8990	112.6735	124.9079	229.5196
18	10.0000	98.5398	112.6735	157.2671	229.5196
19	10.0000	94.5398	112.6735	207.2671	229.5196
20	40.0000	111.9911	112.6735	209.8158	229.5196
21	29.4513	98.5398	112.6735	209.8158	229.5196
22	10.0000	96.7800	108.8846	209.8158	179.5196
23	10.0000	98.5398	68.8846	209.8158	139.7598
24	10.0000	73.4244	30.0000	209.8158	139.7598

Table 3. TABLE III: Summary results for the 10-unit system without loss

Method	Total generation cost ($)			S-time(min)
	Minimum	Average	Maximum
SQP [3]	1051163	NA	NA	0.42
EP [3]	1048638	NA	NA	15.05
CDE [16]	1036756	1040586	1452558	0.20
GA [12]	1033481	1038014	1042606	3.59
EP-SQP [3]	1031746	1035748	NA	7.26
AIS-SQP[28]	1029900	NA	NA	NA
MHEP-SQP[27]	1028924	1031179	NA	21.23
DGPSO [4]	1028835	1030183	NA	4.81
PSO [12]	1027679	1031716	1034340	3.85
SOA [29]	1023946	1026289	1029213	NA
IPSO [8]	1023807	1026863	NA	0.05
CSDE [17]	1023432	1026475	1027634	0.3
CE [14]	1022702	1024024	NA	0.33
ECE [14]	1022272	1023334	NA	0.33
AIS [13]	1021980	1023156	1024973	25.35
ABC[12]	1021576	1022686	1024316	3.47
CDBCO [23]	1021500	1024300	NA	0.73
SOA-SQP[29]	1021460	1023841	1026852	NA
AHDE [10]	1020082	1022474	1024484	1.10
CDE[16]	1019123	1020870	1023115	0.32
HHS [15]	1019091	NA	NA	10.19
ICPSO [9]	1019072	1020027	NA	0.35
CSAPSO[11]	1018767	1019874	NA	0.350
EAPSO [18]	1018510	1018701	1019302	0.63
HIGA [30]	1018473	1019328	1022284	4.41
ICA [20]	1018467	1019291	1021796	NA
TVAC-IPSO[19]	1018217	1018965	1020418	2.72
HBPSO[31]	1018159	1019850	1021813	3.09
CSO [25]	1017660	1018120	1019286	0.90
EBSO [21]	1017147	1017526	1017891	0.15
MILP	1016316			0.94
MILP-IPM	1016311			1.02
NA denotes that the value was not available in the literature.

Table 4. TABLE IV: Outputs (MW) for the 10-unit system without loss

$t$	unit $1$	unit $2$	unit $3$	unit $4$	unit $5$	unit $6$	unit $7$	unit $8$	unit $9$	unit $10$
1	150.0000	142.2665	186.8267	60.0000	122.8666	122.4498	129.5904	47.0000	20.0000	55.0000
2	150.0000	222.2665	230.6932	60.0000	73.0000	122.4498	129.5904	47.0000	20.0000	55.0000
3	150.0000	302.2665	298.6932	60.0000	73.0000	122.4498	129.5904	47.0000	20.0000	55.0000
4	226.6242	316.7994	305.6696	60.0000	122.8666	122.4498	129.5904	47.0000	20.0000	55.0000
5	226.6242	396.7994	299.6696	60.0000	122.8666	122.4498	129.5904	47.0000	20.0000	55.0000
6	303.2484	396.7994	297.3995	60.0000	172.7331	146.2292	129.5904	47.0000	20.0000	55.0000
7	379.8726	396.7994	303.3995	66.8430	172.7331	122.4498	129.5904	55.3121	20.0000	55.0000
8	379.8726	396.7994	297.3995	116.8430	172.7331	122.4498	129.5904	85.3121	20.0000	55.0000
9	456.4968	396.7994	297.3995	131.6073	222.5997	122.4498	129.5904	90.0000	22.0571	55.0000
10	456.4968	396.7994	297.8493	181.6073	222.5997	160.0000	129.5904	120.0000	52.0571	55.0000
11	456.4968	396.7994	321.8493	231.6073	222.5997	160.0000	129.5904	120.0000	52.0571	55.0000
12	456.4968	460.0000	332.6487	231.6073	222.5997	160.0000	129.5904	120.0000	52.0571	55.0000
13	456.4968	396.7994	297.8493	181.6073	222.5997	160.0000	129.5904	120.0000	52.0571	55.0000
14	456.4968	396.7994	297.3995	131.6073	222.5997	122.4498	129.5904	90.0000	22.0571	55.0000
15	379.8726	396.7994	312.5546	110.0000	172.7331	122.4498	129.5904	77.0000	20.0000	55.0000
16	303.2484	396.7994	297.0454	60.0000	122.8665	122.4498	129.5904	47.0000	20.0000	55.0000
17	226.6242	396.7994	299.6696	60.0000	122.8666	122.4498	129.5904	47.0000	20.0000	55.0000
18	303.2484	396.7994	306.0237	66.8430	172.7331	122.4498	129.5904	55.3121	20.0000	55.0000
19	379.8726	396.7994	297.3995	116.8430	172.7331	122.4498	129.5904	85.3121	20.0000	55.0000
20	456.4968	460.0000	332.5857	120.4152	222.5997	160.0000	129.5904	85.3121	50.0000	55.0000
21	456.4968	389.5329	297.3995	119.5092	222.5997	148.5593	129.5904	85.3121	20.0000	55.0000
22	379.8726	309.5329	283.9997	69.5092	172.7331	122.4498	129.5904	85.3121	20.0000	55.0000
23	303.2484	229.5329	203.9997	60.0000	122.8665	122.4498	129.5904	85.3121	20.0000	55.0000
24	226.6242	222.2665	189.7569	60.0000	73.0000	122.4498	129.5904	85.3121	20.0000	55.0000

Table 5. TABLE V: Summary results for the 5-unit system with loss

Method	Total generation cost ($)			S-time(min)
	Minimum	Average	Maximum
SA[5]	47356	NA	NA	4.40
GA [12]	44862	44922	45894	4.43
APSO[6]	44678	NA	NA	NA
AIS [13]	44385	44759	45554	5.33
PSO [12]	44253	45657	46403	4.73
ABC [12]	44046	44065	44219	4.39
EAPSO[18]	43784	NA	NA	0.42
HBPSO[31]	43223	43732	44252	1.54
DE [7]	43213	43813	44247	6.00
HHS [15]	43155	NA	NA	2.33
TVAC-IPSO[19]	43137	43186	43302	1.07
HIGA [30]	43125	43162	43259	2.06
ICA [20]	43117	43144	43210	NA
IPM	43443			0.05
MILP-IPM	43084			0.87

Table 6. TABLE VI: Dispatch results (MW) for the 5-unit system

$t$	unit $1$	unit $2$	unit $3$	unit $4$	unit $5$	Loss	$\Delta P_t$
1	20.6080	98.5398	30.0000	124.9079	139.7598	3.8155	0.0000
2	10.0000	97.8835	66.5747	124.9079	139.7598	4.1259	0.0000
3	10.0000	98.5398	106.5747	124.9079	139.7598	4.7822	0.0000
4	10.0960	98.5398	112.6735	124.9079	189.7598	5.9770	0.0000
5	10.0000	87.5816	112.6735	124.9079	229.5196	6.6825	0.0001
6	40.0000	99.9508	112.6735	133.7152	229.5196	7.8590	0.0001
7	10.0000	98.5398	112.6735	183.7152	229.5196	8.4480	0.0001
8	12.7090	98.5398	112.6735	209.8158	229.5196	9.2577	0.0000
9	42.7090	105.4824	112.6735	209.8158	229.5196	10.2003	0.0000
10	64.0108	98.5398	112.6735	209.8158	229.5196	10.5595	0.0000
11	75.0000	104.0359	112.6735	209.8158	229.5196	11.0448	0.0000
12	75.0000	124.7111	112.6735	209.8158	229.5196	11.7200	0.0000
13	64.0108	98.5398	112.6735	209.8158	229.5196	10.5595	0.0000
14	49.6196	98.5398	112.6735	209.8158	229.5196	10.1683	0.0000
15	19.6196	98.5398	112.6735	202.8589	229.5196	9.2113	0.0001
16	10.0000	82.1494	112.6735	152.8589	229.5196	7.2013	0.0001
17	10.0000	87.5816	112.6735	124.9079	229.5196	6.6825	0.0001
18	10.0000	98.5398	112.6735	165.2180	229.5196	7.9509	0.0000
19	12.7090	98.5398	112.6735	209.8158	229.5196	9.2577	0.0000
20	42.7090	119.9393	112.6735	209.8158	229.5196	10.6572	0.0000
21	39.3529	98.5398	112.6735	209.8158	229.5196	9.9016	0.0000
22	10.0000	98.5398	112.6735	209.8158	181.8844	7.9136	0.0001
23	12.4371	98.5398	72.6735	209.8158	139.7598	6.2261	0.0001
24	10.0000	75.8153	32.6735	209.8158	139.7598	5.0644	0.0000

Table 7. TABLE VII: Summary results for the 10-unit system with loss

Method	Total generation cost ($)			S-time(min)
	Minimum	Average	Maximum
EP[3]	1054685	1057323	NA	47.23
EP-SQP[3]	1052668	1053771	NA	27.53
GA[12]	1052251	1058041	1062511	4.59
MHEP-SQP[27]	1050054	1052349	NA	24.33
DGPSO [4]	1049167	1051725	NA	5.99
PSO[12]	1048410	1052092	1057170	5.45
IPSO [8]	1046275	1048145	NA	0.15
AIS [13]	1045715	1047050	1048431	30.96
CE [14]	1044051	1045159	NA	0.60
ECE [14]	1043989	1044470	NA	0.64
ABC [12]	1043381	1044963	1046805	4.55
CDBCO[23]	1042900	1044700	NA	1.66
HIGA [30]	1041088	1042980	1044927	4.75
TVAC-IPSO[19]	1041066	1042118	1043625	3.16
ICA [20]	1040758	1041665	1043175	NA
EBSO [21]	1038915	1039188	1039272	0.17
CSO[25]	1038320	1039374	1042518	1.39
CSAPSO [11]	1038251	1039543	NA	0.83
EAPSO[18]	1037898	1038109	1038238	2.88
IPM	1047294			0.22
MILP-IPM	1040676			1.12

Table 8. TABLE VIII: Dispatch results (MW) for the 10-unit system

$t$	unit $1$	unit $2$	unit $3$	unit $4$	unit $5$	unit $6$	unit $7$	unit $8$	unit $9$	unit $10$	Loss	$\Delta P_t$	$\Delta P_t$[25]	$\Delta P_t$[18]
1	150.0000	142.2665	199.1034	60.0000	122.8666	122.4498	129.5904	47.0000	20.0000	55.0000	12.2767	0.0000	0.0030	0.0000
2	150.0000	222.2665	238.5251	60.0000	73.0000	122.4498	129.5904	55.3121	20.0000	55.0000	16.1439	0.0000	0.0011	0.0001
3	150.0000	302.2665	312.1763	60.0000	73.0000	122.4498	99.5904	85.3121	20.0000	55.0000	21.7952	0.0001	0.0054	0.0254
4	226.6242	312.1330	297.3995	60.0000	122.8665	122.4498	129.5904	85.3121	20.0000	55.0000	25.3756	0.0001	0.0001	0.0001
5	226.6242	392.1330	297.8235	60.0000	122.8666	122.4498	129.5904	85.3121	20.0000	55.0000	31.7997	0.0001	0.0022	0.0001
6	303.2484	396.7994	297.3995	60.0000	172.7331	144.3175	129.5904	85.3121	20.0000	55.0000	36.4004	0.0000	0.0015	0.0001
7	379.8726	396.7994	306.5561	70.4152	172.7331	122.4498	129.5904	85.3121	22.0571	55.0000	38.7859	0.0001	0.0087	0.1612
8	379.8726	396.7994	302.5307	120.4152	172.7331	122.4498	129.5904	85.3121	52.0571	55.0000	40.7605	0.0001	0.0014	0.0000
9	456.4968	396.7994	297.3995	148.8859	222.5997	122.4498	129.5904	90.0000	52.0571	55.0000	47.2787	0.0001	0.0023	0.0000
10	456.4968	396.7994	337.6226	191.2457	222.5997	160.0000	129.5904	120.0000	52.0571	55.0000	49.4117	0.0000	0.0066	0.0000
11	460.4568	396.7994	340.0000	241.2457	243.0000	160.0000	129.5904	120.0000	52.0571	55.0000	52.1495	0.0001	0.0011	0.0000
12	456.4968	460.0000	340.0000	241.2457	237.0093	160.0000	129.5904	120.0000	80.0000	55.0000	59.3423	0.0001	0.0023	0.0001
13	456.4968	396.7994	337.6226	191.2457	222.5997	160.0000	129.5904	120.0000	52.0571	55.0000	49.4117	0.0000	0.0005	0.0002
14	456.4968	396.7994	297.3995	148.8859	222.5997	122.4498	129.5904	90.0000	52.0571	55.0000	47.2787	0.0001	0.0044	0.0000
15	379.8726	396.7994	314.9562	110.0000	172.7331	122.4498	129.5904	85.3121	50.0000	55.0000	40.7137	0.0001	0.0167	0.0002
16	303.2484	396.7994	292.4895	60.0000	122.8665	122.4498	129.5904	85.3121	20.0000	55.0000	33.7562	0.0001	0.0025	0.0000
17	226.6242	396.7994	324.8902	60.0000	122.8666	122.4498	129.5904	55.3121	20.0000	55.0000	33.5327	0.0000	0.0113	0.0001
18	303.2484	396.7994	306.5303	70.4152	172.7331	122.4498	129.5904	85.3121	22.0571	55.0000	36.1359	0.0001	0.0048	0.0000
19	379.8726	396.7994	302.5307	120.4152	172.7331	122.4498	129.5904	85.3121	52.0571	55.0000	40.7605	0.0001	0.0047	0.0002
20	456.4968	460.0000	340.0000	135.3837	222.5997	160.0000	129.5904	115.3121	52.0571	55.0000	54.4398	0.0000	0.0038	0.0000
21	456.4968	389.5329	302.4598	120.4152	222.5997	160.0000	129.5904	85.3121	50.0000	55.0000	47.4070	0.0001	0.0060	0.0000
22	379.8726	309.5329	297.3995	87.6310	172.7331	122.4498	129.5904	85.3121	20.0000	55.0000	31.5216	0.0002	0.0070	0.0000
23	303.2484	229.5329	224.1642	60.0000	122.8666	122.4498	129.5904	85.3121	20.0000	55.0000	20.1645	0.0001	0.0075	0.0001
24	226.6242	222.2665	205.8519	60.0000	73.0000	122.4498	129.5904	85.3121	20.0000	55.0000	16.0950	0.0001	0.0001	0.0002

Table 9. TABLE IX: Solution differences between MILP and MILP-IP

case 3: 5-unit system
0	(0,3]	(3,6]	(6,9]	(7,12]
66.67%	12.50%	8.33%	7.50%	5.00%
case 4: 10-unit system
0	(0,10]	(10,20]	(20,30]	(30,40]
70.83%	15.00%	5.42%	5.00%	3.75%

Equations

$$\min \sum _{t=1}^T\sum _{i=1}^{N}c(P_{i,t})$$

$$c(P_{i,t})={\alpha }_i+{\beta }_i{P}_{i,t}+{\gamma }_i{P}_{i,t}^2+e_i|\sin (f_i(P_{i,t}-P_i^{\min }))|$$

$$\sum _{i=1}^N{P}_{i,t}=D_t+P_t^{loss},\forall t$$

$$P_t^{loss}=\sum _{i=1}^N\sum _{j=1}^N{P}_{i,t}{B}_{i,j}{P}_{j,t},\forall t$$

$$P_i^{\min }\le P_{i,t}\le P_i^{\max },\forall i,t$$

$$D{R}_i\le P_{i,t}-P_{i,t-1}\le U{R}_i,\forall i,t$$

$$\{\begin{array}{rl}& S{R}_{i,t}\le P_i^{\max }-P_{i,t},\\ & S{R}_{i,t}\le \tau U{R}_i,\forall i,t\end{array}$$

$$\sum _{i=1}^{N}S{R}_{i,t}\ge R_t,\forall t$$

$$\min$$

$$s.t.$$

$$s_{i,t}\ge -\sin (f_i(P_{i,t}-P_i^{\min }))$$

$$s_{i,t}-\sin (f_i(P_{i,t}-P_i^{\min }))-u_{i,t}^0=0$$

$$s_{i,t}+\sin (f_i(P_{i,t}-P_i^{\min }))-v_{i,t}^0=0$$

$$2s_{i,t}-u_{i,t}^0-v_{i,t}^0=0$$

$$-2\sin (f_i(P_{i,t}-P_i^{\min }))-u_{i,t}^0+v_{i,t}^0=0$$

$$s_{i,t}-\frac{u_{i,t}^0}{2}-\frac{v_{i,t}^0}{2}=0$$

$$\sin (f_i(P_{i,t}-P_i^{\min }))+\frac{u_{i,t}^0}{2}-\frac{v_{i,t}^0}{2}=0$$

$$s_{i,t}-u_{i,t}-v_{i,t}=0$$

$$\sin (f_i(P_{i,t}-P_i^{\min }))+u_{i,t}-v_{i,t}=0$$

$$u_{i,t}\ge 0,v_{i,t}\ge 0$$

$$\min$$

$$s.t.$$

$$\min$$

$$s.t.$$

$$a_{l-1,i}{z}_{l,i,t}\le P_{l,i,t}\le a_{l,i}{z}_{l,i,t}$$

$$\sum _{l=1}^{L_i}{z}_{l,i,t}=1$$

$$z_{l,i,t}\in \{0,1\}$$

$$(\text{ref}balance),(\text{ref}G-limits)-(\text{ref}SR2)$$

$$\{\begin{array}{rl}& L_i=ceil(M\frac{f_i(P_i^{max}-P_i^{min})}{\pi })\\ & k_{l,i}=\frac{c(a_{l,i})-c(a_{l-1,i})}{a_{l,i}-a_{l-1,i}}\\ & b_{l,i}=c(a_{l-1,i})-k_{l,i}{a}_{l-1,i}.\end{array}$$

$$ScaledCPUtime$$

$$=\frac{GivenCPUspeed}{BaseCPUspeed}GivenCPUtime.$$

$$\Delta P_t=|\sum _{i=1}^N{P}_{i,t}-D_t-\sum _{i=1}^N\sum _{j=1}^N{P}_{i,t}{B}_{i,j}{P}_{j,t}|.$$

Full text

I Introduction

Dynamic economic dispatch (DED) is one of the fundamental issues for the optimal economic operation in power system, which aims at allocating the customers’ load demands among the available thermal power generating units in an economic, secure and reliable way at a certain time of interest [1]. Traditionally, the generation cost function of DED is approximated by a convex, quadratic and differentiable polynomials. However, in actual operation, wire drawing effects, occurring as each steam admission valve in a turbine starts to open, produce a rippling effect on the generation cost curve[2], which is known as the valve-point effect (VPE). To model the effect of valve-points, a recurring rectified sinusoid contribution is added to the input-output equation [2], which makes the generation cost function non-convex, non-differentiable and multiple extremal. When VPE is ignored, the rough approximation of the generation cost function will introduce some inaccuracies into the dispatch results. In order to improve the optimality of the solution, a more accurate DED model, dynamic economic dispatch with valve-point effect (DED-VPE) should be considered. However, when VPE is considered, some non-convex, non-differentiable and multiple extremal characteristics are introduced, which makes the solution for DED-VPE more challenging.

In the past decades, a number of optimization methods have been proposed to solve the DED-VPE. Due to the intractability of the problem, most of the currently available approaches for DED-VPE are heuristic optimization techniques[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26], such as genetic algorithm (GA) [2], evolutionary programming (EP) [3], simulated annealing (SA) [5], particle swarm optimization (PSO) [6], differential evolution algorithm (DE) [7], artificial bee colony algorithm (ABC) [12], artificial immune system (AIS) [13], enhanced cross-entropy (ECE) [14], harmony search (HS) [15], imperialist competitive algorithm (ICA) [20], bee swarm optimization algorithm (BSO) [21], teaching-learning algorithm (TLA) [24], etc. These heuristic optimization techniques are population-based search methods which do not depend on the gradient and Hessian operators of the objective function. So, they can be applied to solve the DED-VPE problem effectively. However, they are quite sensitive to various parameter settings and solution may be different at each trial. Hence, hybrid methods which combine several heuristic techniques or deterministic approaches [3, 27, 28, 29, 30, 31] such as hybrid evolutionary programming and sequential quadratic programming (EP-SQP) [3], hybridization of artificial immune systems and sequential quadratic programming (AIS-SQP) [28], hybrid seeker optimization algorithm and sequential quadratic programming (SOA-SQP) [29], hybrid immune-genetic algorithm (HIGA) [30], etc, tend to be more efficient than the individual methods. However, they still have the intrinsic drawbacks of the heuristic method we mentioned above.

Unlike heuristics, deterministic mathematical programming-based optimization techniques can solve to a robust result due to the solid mathematical foundations and the availability of the powerful software tools. Therefore, a strategy recently appeared for DED-VPE is to reformulate the generation cost function, yielding a good optimization model that can be solved by a deterministic method. In [32], the non-convex and non-differentiable cost function caused by VPE is piecewise linearized, then a mixed integer quadratic programming (MIQP) method can be used to solve DED-VPE. But when the MIQP formulation is directly solved by using a mixed integer programming (MIP) solver, the optimization will suffer convergence stagnancy and run out of memory. As a result, the multi-step method, the warm start technique and the range restriction scheme are required. However, the range restriction scheme just restricts the solution space to a subspace where the global optimal solution would probably lie in. Consequently, the optimality of the solution for the MIQP can not be guaranteed. Besides, when the complicated transmission loss is considered, base on the above process, more tedious adjustment techniques are needed. Different from [32], the whole generation cost function is considered for piecewise linearization in [33], then a mixed integer linear programming (MILP) formulation is proposed to solve the DED-VPE, where a global optimal solution within a preset tolerance can be guaranteed. But the transmission loss is not considered in [33]. Therefore, more efforts are worth developing effective deterministic mathematical programming-based optimization methods for DED-VPE to obtain better dispatch results.

In this paper, a hybrid deterministic method that integrates the MILP and interior point method (IPM), denoted by MILP-IPM, is proposed for solving the DED-VPE problem while transmission loss is included. Due to the non-differentiable characteristic of DED-VPE, the classical derivative-based optimization methods can not be used any more. With the help of model reformulation, we derive a non-linear programming (NLP) formulation of DED-VPE, which can be solved by the polynomial time IPM immediately. However, IPM is a local optimization method. If the DED-VPE is solved by IPM in a single step, the optimization will easily trap in a poor local optima due to its non-convex and multiple local minima characteristics. In order to overcome this deficiency, MILP method [33] is combined to generate a good initial point for IPM. And then, solving its NLP formulation via IPM, a good optimal solution for DED-VPE can thus be determined.

The rest of this paper is organized as follows. Section 2 describes the mathematical formulation of DED-VPE. Section 3 derives an NLP formulation for DED-VPE. Section 4 introduces the implementation of MILP-IPM. Section 5 presents the simulation results and analysis. Section 6 draws the conclusions.

II Mathematical formulation of DED-VPE

The DED-VPE problem usually can be formulated as a non-convex and non-differentiable optimization problem. The objective of DED-VPE is to minimize the total generation cost over a scheduled time horizon, which can be written as:

[TABLE]

where

[TABLE]

where $P_{i,t}$ is the power output of uint $i$ in period $t$; $P_{i}^{\min}$ is the minimum power output of unit $i$; $\alpha_{i}$, $\beta_{i}$, $\gamma_{i}$, $e_{i}$ and $f_{i}$ are positive coefficients of unit $i$; $N$ and $T$ are total number of units and periods, respectively.

The minimized DED-VPE problem should be subjected to the constraints as follows.

• •

Power balance equations

[TABLE]

where $D_{t}$ is the load demand in period $t$; $P_{t}^{loss}$ is the transmission loss in period $t$, which can be calculated based on the B-Matrix loss coefficients and expressed in the quadratic form as given below [34]

[TABLE]

where $B_{i,j}$ is the $(i,j)$-th element of the transmission loss coefficient matrix.

• •

Power generation limits

[TABLE]

where $P_{i}^{\max}$ is the maximum power output of unit $i$.

• •

Ramp rate limits

[TABLE]

where $DR_{i}$ and $UR_{i}$ are the ramp-down and ramp-up rates of unit $i$, respectively.

• •

Spinning reserve constraints

[TABLE]

[TABLE]

where $SR_{i,t}$ is the spinning reserve provided by unit $i$ in period $t$; $R_{t}$ is the system spinning reserve requirement in period $t$; $\tau$ is the time duration for units to deliver reserve [22, 32].

III An NLP formulation for DED-VPE

As we can see from the section II, DED-VPE is a non-convex and non-differentiable optimization problem which is hard to tackle. Due to the non-differentiable characteristic of DED-VPE, the classical mathematical programming-based methods, also known as derivative-based optimization methods, are not suitable any more. To overcome this difficulty, we replace $|\sin(f_{i}(P_{i,t}-P_{i}^{\min}))|$ for (2) with an auxiliary variable $s_{i,t}$, then the objective function (1) can be equivalent to

[TABLE]

By introducing some slack variables $u_{i,t}^{0}\geq 0$, $v_{i,t}^{0}\geq 0$, the new inequality constraints (10) and (11) are converted into equality constraints

[TABLE]

From (12) and (13), we can get the following equations

[TABLE]

Then we have

[TABLE]

Let $u_{i,t}=\frac{u_{i,t}^{0}}{2}$, $v_{i,t}=\frac{v_{i,t}^{0}}{2}$, we can obtain that

[TABLE]

Consequently, the DED-VPE can be formulated as the following differentiable NLP formulation.

[TABLE]

IV Implementation of MILP-IPM

It is well known that, IPM is a powerful tool for solving non-linear optimization problem and it has successfully been implemented in various issues in power system [35], such as optimal power flows, state estimation, hydrothermal coordination, economic dispatch, etc. Hence, the above NLP formulation (21) for DED-VPE can be directly solved by using IPM. However, if the NLP formulation of DED-VPE is solved by IPM in a single step, the solution will easily trap in a poor local optima due to its non-convex, non-linear and multiple local minima characteristics.

To exploit a better solution, an MILP method presented in [33] is incorporated to find a good initial point for IPM. In [33], when transmission loss is not considered, the DED-VPE is reformulated into an MILP formulation which can be solved by a state-of-the-art MIP solver directly and efficiently. Hence, a global optimal solution within a preset tolerance, which can be used as an initial point for IPM to improve the quality of the eventual dispatch result, can be guaranteed via an enumeration algorithm.

IV-A An MILP formulation for DED-VPE

To obtain an MILP formulation of DED-VPE, ${L_{i}}+1$ break points are chosen over a generation interval $[P_{i}^{\min},P_{i}^{\max}]$, such that $P_{i}^{min}=a_{0,i}\leq a_{1,i}\leq\cdots\leq a_{{L_{i}},i}=P_{i}^{max}$. Segment variables $P_{l,i,t}$ and binary variables $z_{l,i,t}$ $(l=1,\cdots{L_{i}})$ are introduced to make $P_{m,i,t}=P_{i,t}$ and $P_{l,i,t}=0(l\neq m)$ when the $P_{i,t}$ lies in segment $m(m\in\{1,\cdots{L_{i}}\})$. The detail of the MILP formulation for DED-VPE without considering the transmission loss is expressed as [33]:

[TABLE]

where $P_{t}^{loss}=0$ for (3) and $L_{i}$, $k_{l,i}$, $b_{l,i}$ are calculated as follows

[TABLE]

Above $ceil(x)$ means round $x$ to the nearest integer greater than or equal to $x$ and $M$ is the number of equal segments on each $sin(x)$ where $x$ belongs to $[0,\pi]$.

IV-B Details of MILP-IPM

Now the MILP-IPM for DED-VPE as a whole is summarized as follows.

Step 1: Solve the MILP formulation (22) by using MILP method to obtain a global optimal solution within a preset tolerance for DED-VPE without transmission loss.

Step 2: Solve the NLP formulation (21) by using IPM, where the initial point is taken at the solution gained in Step 1, to obtain a good local optimal solution for DED-VPE with transmission loss.

In step 1, an MILP formulation for DED-VPE without transmission loss is solved, yielding a solution which is close to a good optimal solution for DED-VPE with transmission loss. This is mainly because, in a DED problem, the transmission loss at each period is small compared to the corresponding load demand [32]. When transmission loss is ignored, the solution gained in step 1 is “the most” economic. When transmission loss is considered, it means that more outputs are needed and some transmission loss constraints will be added into the original model. Then base on the most economic initial solution, some units outputs will be fine tuned (since the transmission loss is small compared to the load demand ) by IPM in step 2, to meet the new constraints and attain a new economic solution. The following simulation results conform well with the above analysis.

V Simulation results and analysis

To assess the validity and effectiveness of the proposed MILP-IPM, in this section, several test systems which are widely studied for DED-VPE with and without consideration of transmission loss constraints over a scheduled time horizon of 24 h are simulated. For a fairer comparison with most of the existing methods, the spinning reserve constraints (7) and (8) are not included in our simulations. Actually, it does not reduce the difficulty of the problem since the spinning reserve constraints are some simple linear constraints which is very easy to handle for MILP and IPM.

Since the computation time highly depends on the computer system used, so the CPU execution times for different methods may not be directly comparable due to different computers used. In order to have a fair comparison regarding the computational effort, the CPU chip frequency from the used computer is used to convert the CPU times obtained from different methods into a common base for comparison [36]:

[TABLE]

In this paper, the base CPU speed is 2.4 GHz and we denote $Scaled~{}CPU~{}time$ by S-time for short.

In this section, we first carry out two experiments on DED-VPE without transmission loss to compare MILP-IPM with MILP and other methods, showing the potential of MILP-IPM for seeking better solutions. Then two experiments on DED-VPE with transmission loss are carried out to demonstrate the validity and effectiveness of MILP-IPM in solving DED-VPE. All cases are performed on an Intel Core 2.5 GHz Dell-notebook with 8 GB of RAM. Our models are coded in YALMIP [37] within the MATLAB R2014a and optimised using CPLEX 12.6.1 [38] for solving the MILP (22) and IPOPT 3.12.6 [39] for solving the NLP (21).

V-A DED-VPE without transmission loss

Case 1: A 5-unit system without transmission loss

The first test case is incorporated five thermal units. The characteristics of the thermal units and load demands are taken from [5]. Firstly, we directly solve its MILP formulation using CPLEX to 3.2% optimality. In this formulation, segment parameter $M$ is set to 4. After the optimization, an optimal solution is yielded in 0.82 min with a total generation cost 42563\$$. Taking such a solution as an initial point in step 2, we solve its NLP formulation using IPOPT with the default options, then a local optimal solution with an optimal value 42524$$ is found in 0.04 min.

Table I lists comparison results of the total generation cost obtained by MILP and MILP-IPM. It is clearly seen that the proposed MILP-IPM can make a certain improvement for 5-unit system.

The outputs obtained by using MILP-IPM for the 5-unit system are given in Table II for verification.

Case 2: A 10-unit system without transmission loss

The second test case contains ten thermal units. The characteristics of the thermal units and load demands are taken from [3]. Firstly, we directly solve its MILP formulation using CPLEX to 0.3% optimality. After the optimization, an optimal solution is yielded in 0.94 min with a total generation cost 1016316\$$. Using such a solution as an initial point in step 2, we solve its NLP formulation using IPOPT with the default options, then a local optimal solution with optimal value 1016311$$ is found in 0.08 min.

The comparison results of the total generation cost obtained by using MILP-IPM and other methods are shown in Table III. It is obvious that the proposed MILP-IPM can solve to the lowest generation cost among all the methods in a reasonable time. But at the same time, we also see that, although MILP-IPM can exploit a better solution in comparison with MILP, but the improvement is small. This is mainly because the solution achieved in step 1 is very well.

The outputs obtained by using MILP-IPM for the 10-unit system are given in Table IV for verification.

V-B DED-VPE with transmission loss

Since transmission loss can not be avoided in a power distribution system and it is critical in real-world DED problem, solution for DED-VPE with transmission loss has more values. Due to the data unavailability [26], in our simulations, only two cases (5- and 10- unit systems) for DED-VPE with transmission loss are considered.

Case 3: A 5-unit system with transmission loss

The third test case is a 5-unit system in which the characteristics of the thermal units and load demands are the same as those in case 1. In this case, the transmission loss is considered. Owing to the limits of space, the loss coefficients are not listed here. One can refer to [5].

Since the transmission loss is included, unlike case 1, it can not be solved by MILP immediately. Fortunately, with the help of model reformulation, we derive a differentiable NLP formulation (21) of DED-VPE, which can be solved by the powerful IPM directly.

In order to demonstrate the validity and effectiveness of the proposed NLP formulation and MILP-IPM for the DED-VPE with transmission loss, on the one hand, we directly solve its NLP formulation using IPOPT with the default options. Then a local optimal solution with an optimal value 43443\$$ is obtained in 0.05 min. On the other hand, we solve this system by using MILP-IPM, where the parameters in step 1 are set the same as the case 1. After the optimization, a local optimal solution with an optimal value 43084$$ is found in 0.87 min. The results are compared with other methods in Table V.

As we can see in Table V, although the proposed NLP formulation can be solved to a local optimal solution in a short time, but when the DED-VPE is directly solved by IPM in a single step, the obtaining solution is not the best due to its non-convex and multiple local minima characteristics. While the MILP-IPM which combines MILP with IPM is employed, it can solve to the lowest generation cost among all the methods in a reasonable time.

But, we should note that, not only the optimality but also the feasibility should be considered for assessing the quality of the solution. In Table VI, the generation dispatch results obtained by using MILP-IPM for the 5-unit system are presented, where $\Delta P_{t}$ denotes the violation degree of the power balance constraint at interval $t$, which is calculated by

[TABLE]

Note that, the outputs which have been round-off are used for this calculation (i.e. the outputs in Table VI). Actually, when the original outputs are adopted, all the violations are less than 7e-7. In other words, the solution obtained by using MILP-IPM strictly satisfies the power balance constraints at the same time. This is mainly the result of the rigorous theoretical foundations of interior point algorithm.

Case 4: A 10-unit system with transmission loss

The fourth test case is a 10-unit system in which the characteristics of the thermal units and load demands are the same as those in case 2 and the loss coefficients are taken from [4].

Similar to the case 3, on the one hand, we directly solve its NLP formulation using IPOPT with the default options. Then a local optimal solution with an optimal value 1047294\$$ is found in 0.22 min. On the other hand, we solve this system by using MILP-IPM, where the parameters in step 1 are set the same as the case 2. After the optimization, a local optimal solution with an optimal value 1040676$$ is found in 1.12 min. The results are compared with other methods in Table VII.

As we can see in the Table VII, in comparison with IPM, MILP-IPM can solve to a much better solution in a reasonable time. Meanwhile, the total generation cost obtained by MILP-IPM is lower than most of the results reported in the literatures.

In Table VIII, we provide the generation dispatch results obtained by MILP-IPM for the 10-unit system. In this table, the violations of power balance constraints for CSO[25] and EAPSO[18] are also calculated due to their available solutions. For a fair comparison, we use the outputs shown in the corresponding literatures to compute the corresponding $\Delta P_{t}$. From the Table VIII we notice that, although CSO[25] and EAPSO[18] can obtain lower total generation costs than MILP-IPM, but MILP-IPM outperforms them in terms of the feasibility of the solution for DED-VPE. In fact, when the original outputs are adopted for our case, all the violations are less than 8e-6. It means that, the feasibility of the solution obtained by MILP-IPM can be strictly satisfied.

V-C Results analysis

From the simulation process in this section, we observe that, in the 5-unit system with transmission loss, $66.67\%$ of the solution points obtained by MILP-IPM are the same as those obtained by the individual MILP and differences of the rest range from 0.0960 to 11.72. In the 10-unit system with transmission loss, $70.83\%$ of the solution points obtained by MILP-IPM are the same as those obtained by the individual MILP and differences of the rest range from 0.0001 to 39.77. The details with respect to the solution differences between MILP and MILP-IP for both cases are given in Table IX.

To some extent, this phenomenon indicates that after solving MILP formulation in step 1, an initial point which is close to a good optimal solution for DED-VPE is yielded, and then in step 2, IPM starts its search from this initial point and tunes to a good optimal solution where all the constraints are satisfied, which is consistent with the analysis in the subsection 4.2.

VI Conclusion

In this paper, a deterministic MILP-IPM is proposed to solve the non-convex and non-differentiable DED-VPE. To avoid the intractable non-differentiable characteristic of DED-VPE, we derive a differentiable NLP formulation for DED-VPE. Although the NLP formulation can be directly solved by IPM, but the optimization will easily trap in a poor local optima. Therefore, MILP is integrated to generate a good initial point. And then, IPM can be used to exploit a better local optima for DED-VPE. Comparing with the heuristics which are inherently stochastic, MILP-IPM results are much more stable and the feasibility of the solution can be strictly guaranteed. So, MILP-IPM as a deterministic optimization technique is very promising to apply to the practical problems when VPE type factors are considered.