Pareto Optimal Demand Response Based on Energy Costs and Load Factor in Smart Grid

TL;DR

This paper introduces a multiobjective optimization approach for demand response in smart grids, balancing residential energy costs and load factors to improve grid efficiency and user satisfaction.

Contribution

It develops a novel multiobjective evolutionary algorithm that effectively generates Pareto optimal demand response solutions considering energy costs and load factors.

Findings

Outperforms comparable methods in reducing energy costs.

Achieves a balanced load factor for the grid.

Produces Pareto optimal demand response programs.

Abstract

Demand response for residential users is essential to the realization of modern smart grids. This paper proposes a multiobjective approach to designing a demand response program that considers the energy costs of residential users and the load factor of the underlying grid. A multiobjective optimization problem (MOP) is formulated and Pareto optimality is adopted. Stochastic search methods of generating feasible values for decision variables are proposed. Theoretical analysis is performed to show that the proposed methods can effectively generate and preserve feasible points during the solution process, which comparable methods can hardly achieve. A multiobjective evolutionary algorithm is developed to solve the MOP, producing a Pareto optimal demand response (PODR) program. Simulations reveal that the proposed method outperforms the comparable methods in terms of energy costs while producing a satisfying load factor. The proposed PODR program is able to systematically balance the needs of the grid and residential users.

Tables

Table 1. TABLE I: Notations

$h$	Time slot $h$ where $h\in \mathscr{H}$
$H$	Observation horizon (the size of $\mathscr{H}$)
$𝒜$	Set of unshiftable and inflexible appliances
$a$	Unshiftable and inflexible appliance
$p_a^h$	Amount of power for appliance $a$ to operate in time slot $h$
${𝒜}_S$	Set of shiftable and inflexible appliances
$b$	Shiftable and inflexible appliance
$h_b$	Working time slot of appliance $b$
$s_b$ and $e_b$	Start and end time slots of appliance $b$
$w_b$	Total working hours of appliance $b$
$p_b^h\left(h_b\right)$	Amount of power for appliance $b$ to operate in time slot $h$
${𝒜}_{SF}$	Set of shiftable and flexible appliances
$c$	Shiftable and flexible appliance
$s_c$ and $e_c$	Start and end time slots of appliance $c$
$p_c^h$	Amount of power for appliance $c$ to operate in time slot $h$
$p_c^{\max }$	Maximum operating power of appliance $c$
$p_c^{\min }$	Minimum operating power of appliance $c$
$p_{c,\text{agr}}^{\min }$	Aggregate minimum power for using appliance $c$
$s_{\text{ev}}$ and $e_{\text{ev}}$	Start and end time slots of EV battery charging
$p_{\text{ev}}^h$	EV charging power in time slot $h$
$p_{\text{ev}}^{\max }$	Maximum charging rate of an EV
$B_{\text{ev}}^0$	Initial capacity of the EV battery
$B_{\text{ev}}^{\min }$	Minimum capacity of the EV battery
$B_{\text{ev}}^{\max }$	Maximum capacity of the EV battery
$B^0$	Initial capacity of the residential energy storage system
$B^h$	Energy level of the residential energy storage system in time slot $h$
$B^{\min }$	Minimum capacity of the residential energy storage system
$B^{\max }$	Maximum capacity of the residential energy storage system
$u^h$	Charging and discharging control for the residential energy storage system in time slot $h$
$r^h$	Amount of power provided by renewable energy sources in time slot $h$

Table 2. TABLE II: Simulation Parameters

Notations	Power	Working Schedule
$p_{a_1}^h$	0.02 kW	$h\in \{17,18,\dots ,24\}$
$p_{a_2}^h$	0.22 kW	unif(18, 22, 3)
$p_{a_3}^h$	0.2 kW	unif(11, 13, 3)
$p_{a_4}^h$	0.2 kW	unif(16, 18, 5)
$p_{a_5}^h$	0.7 kW	unif(18, 22, 1)
$p_{a_6}^h$	1.3 kW	unif(14, 16, 1)
$p_{a_7}^h$	0.2 kW	unif(18, 22, 1)
$p_{a_8}^h$	0.08 kW	unif(18, 20, 3)
$p_{a_9}^h$	0.05 kW	$h\in \{1,2,\dots ,24\}$
$p_{a_{10}}^h$	1.5 kW	$h=8$
$p_{a_{11}}^h$	1.6 kW	$h\in \{17,18\}$
$p_{a_{12}}^h$	0.2 kW	$h\in \{1,2,\dots ,24\}$
$p_{a_{13}}^h$	0.8 kW	$h=17$
Notations	Values
$s_{b_1},e_{b_1}$	Random number from $\{10,11,12,13\}$, $e_{b_1}=s_{b_1}+7$
$p_{b_1}^h,w_{b_1}$	1 kW if $h_{b_1}\cap \left\{h\right\}\ne \varnothing ,w_{b_1}=1$
$s_{b_2},e_{b_2}$	Random number from $\{12,13,14,15\}$, $e_{b_2}=s_{b_2}+4$
$p_{b_2}^h,w_{b_2}$	1 kW if $h_{b_2}\cap \left\{h\right\}\ne \varnothing ,w_{b_2}=2$
$s_{b_3},e_{b_3}$	Random number from $\{13,14,15,16\}$, $e_{b_3}=s_{b_3}+7$
$p_{b_3}^h,w_{b_3}$	2 kW if $h_{b_3}\cap \left\{h\right\}\ne \varnothing ,w_{b_3}=2$
$s_{c_1}$	12
$e_{c_1}$	24
$s_{c_2}$	Random number from $\{20,21,22,23\}$
$e_{c_2}$	$s_{c_2}+9$
$s_{\text{ev}}$	Random number from $\{18,19,\dots ,22\}$
$e_{\text{ev}}$	$s_{\text{ev}}+11$
$p_{c_1}^{\max },p_{c_2}^{\max }$	3 kW
$p_{c_1}^{\min },p_{c_2}^{\min }$	0.5 kW
$p_{c_1,\text{agr}}^{\min }$	29 kW
$p_{c_2,\text{agr}}^{\min }$	12 kW
$p_{\text{ev}}^{\max }$	3 kW
$B^{\max },B^0$	4 kWh, 1 kWh
$B_{\text{ev}}^{\max },B_{\text{ev}}^{\min }$	24 kWh, $B_{\text{ev}}^{\min }=0.8B_{\text{ev}}^{\max }$
$B_{\text{ev}}^0$	Random number from $[0.3B_{\text{ev}}^{\max },0.6B_{\text{ev}}^{\max }]$

Table 3. TABLE III: Comparison of Various Power Scheduling Methods

Cost Increase With Respect to
PODR Program
Date	Load Variance	Load Factor	Area-load	Payment Minimization
	Minimization	Maximization	Method
July 27, 2017	26.4%	30.6%	14.4%	8.5%
July 6, 2017	18.6%	26%	15.5%	2.9%
July 23, 2017	15%	23.1%	14.1%	7.3%
July 29, 2017	12.4%	19.6%	14%	6.7%
average	18.1%	24.8%	14.5%	6.4%
Load Factor Reduction With Respect to
Modified Load Variance Minimization Method
Date	PODR	Load Factor	Area-load	Payment Minimization
	Program	Maximization	Method
July 27, 2017	-13.1%	-4.1%	-10.6%	-24.5%
July 6, 2017	-8.2%	-3.8%	-6.5%	-23.6%
July 23, 2017	-12.7%	-4.4%	-9.3%	-31.4%
July 29, 2017	-8.4%	-4.4%	-4.7%	-18.8%
average	-10.6%	-4.2%	-7.8%	-24.6%

Equations

$$h_b\in C_{w_b}^{s_b\to e_b}\forall b\in {\mathcal{A}}_S.$$

$$C_{w_b}^{s_b\to e_b}=\left\{\left\{15,16\right\},\left\{15,17\right\},\left\{16,17\right\}\right\}$$

$$p_b^h(h_b)$$

$$p_c^{\min }\le p_c^h\le p_c^{\max }\text{mbox}and\sum _{h=s_c}^{e_c}{p}_c^h\ge p_{c,\text{agr}}^{\min }$$

$$(e_c-s_c+1)p_c^{\max }>p_{c,\text{agr}}^{\min }.$$

$$0\le p_{\text{ev}}^h\le p_{\text{ev}}^{\max }\text{mbox}and{B}_{\text{ev}}^{\min }\le B_{\text{ev}}^0+\sum _{h=s_{\text{ev}}}^{e_{\text{ev}}}{p}_{\text{ev}}^h\Delta h\le B_{\text{ev}}^{\max }$$

$$B_{\text{ev}}^0+(e_{\text{ev}}-s_{\text{ev}}+1)p_{\text{ev}}^{\max }\Delta h>B_{\text{ev}}^{\min }.$$

$$B^h=B^{h-1}+(r^h-u^h)\Delta h.$$

$$0\le B^h\le B^{\max }\forall h\in \mathcal{H}$$

$$E_{\text{total}}^h=\max \{\sum _{a\in \mathcal{A}}{p}_a^h+\sum _{b\in {\mathcal{A}}_S}{p}_b^h(h_b)$$

$$+\sum _{c\in {\mathcal{A}}_{SF}}{p}_c^h+p_{\text{ev}}^h-u^h,0\}\Delta h.$$

$$\sum _{h\in \mathcal{H}}{E}_{\text{total}}^h{\lambda }^h.$$

$$\frac{\sum _{h\in \mathcal{H}}{E}_{\text{total}}^h/H}{\underset{h\in \mathcal{H}}{\max }{E}_{\text{total}}^h}$$

$$\begin{array}{rl}\underset{\mathbf{x}}{\min }& f_{\text{cost}}(\mathbf{x})=\sum _{h\in \mathcal{H}}{E}_{\text{total}}^h{\lambda }^h\\ \underset{\mathbf{x}}{\max }& f_{\text{LF}}(\mathbf{x})=\frac{\sum _{h\in \mathcal{H}}{E}_{\text{total}}^h/H}{\underset{h\in \mathcal{H}}{\max }{E}_{\text{total}}^h}\\ \text{subject to }& (\text{ref}e{q}_{c}s{t}_c),(\text{ref}e{q}_{c}s{t}_{e}v),\text{ and }(\text{ref}e{q}_{c}s{t}_{s}torage)\forall c\in {\mathcal{A}}_{SF}\end{array}$$

$$[\phi(h_{b})]_{j}=\left\{\begin{array}[]{ll}1,&\hbox{ if }j+s_{b}-1\in h_{b}\\ 0,&\hbox{otherwise}\end{array}\right.$$

$$x_{\text{new}}^h=\delta x^h+(1-\delta ){x^{\prime }}^h$$

$$\begin{array}{rl}& p_c^h\in [p_c^{\min },p_c^{\max }],{p^{\prime }}_c^h\in [p_c^{\min },p_c^{\max }],\sum _{h=s_c}^{e_c}{p}_c^h\ge p_{c,\text{agr}}^{\min },\text{mbox}and\\ & \sum _{h=s_c}^{e_c}{p^{\prime }}_c^h\ge p_{c,\text{agr}}^{\min }.\end{array}$$

$$\begin{array}{rl}& \delta p_c^h+(1-\delta ){p^{\prime }}_c^h\in [p_c^{\min },p_c^{\max }]\text{mbox}and\\ & \sum _{h=s_c}^{e_c}\delta p_c^h+(1-\delta ){p^{\prime }}_c^h\ge p_{c,\text{agr}}^{\min }\end{array}$$

$$\begin{array}{rl}& B^{h-1}+(r^h-u^h)\Delta h\in [0,B^{\max }]\text{mbox}and\\ & {B^{\prime }}^{h-1}+(r^h-{u^{\prime }}^h)\Delta h\in [0,B^{\max }].\end{array}$$

$$B_{\text{new}}^h=\delta B^h+(1-\delta ){B^{\prime }}^h\text{mbox}and{u}_{\text{new}}^h=\delta u^h+(1-\delta ){u^{\prime }}^h.$$

$$B_{\text{new}}^{h-1}+(r^h-u_{\text{new}}^h)\Delta h\in [0,B^{\max }].$$

$$\begin{array}{rl}& B_{\text{new}}^h\\ =& \delta (B^{h-1}+(r^h-u^h)\Delta h)+(1-\delta )({B^{\prime }}^{h-1}+(r^h-{u^{\prime }}^h)\Delta h)\\ =& B_{\text{new}}^{h-1}+(r^h-u_{\text{new}}^h)\Delta h.\end{array}$$

$$p_c^h:=p_c^h+(p_c^{\max }-p_c^h){\delta }^h$$

$$p_{\text{ev}}^h:=p_{\text{ev}}^h+(p_{\text{ev}}^{\max }-p_{\text{ev}}^h){\delta }^h$$

$$p_{\text{ev}}^h:=\frac{B_{\text{ev}}^{\max }-B_{\text{ev}}^0}{{\sum }_{h=s_{\text{ev}}}^{e_{\text{ev}}}{p}_{\text{ev}}^h\Delta h}{p}_{\text{ev}}^h$$

$$u^h\in [\frac{B^{h-1}+r^h\Delta h-B^{\max }}{\Delta h},\frac{B^{h-1}+r^h\Delta h}{\Delta h}].$$

$$B^h=B^{h-1}+(r^h-u^h)\Delta h.$$

$$\begin{array}{rl}{\mathbf{x}}^{\ast }=& \arg \underset{\mathbf{x}\in X(t_{\max })}{\min }\frac{f_{\text{cost}}(\mathbf{x})-\underset{{\mathbf{x}}^{\prime }\in X(t_{\max })}{\min }{f}_{\text{cost}}({\mathbf{x}}^{\prime })}{L_{\text{cost}}}\\ & +\frac{\underset{{\mathbf{x}}^{\prime }\in X(t_{\max })}{\max }{f}_{\text{LF}}({\mathbf{x}}^{\prime })-f_{\text{LF}}(\mathbf{x})}{L_{\text{LF}}}\end{array}$$

$$\begin{array}{rl}{L}_{\text{cost}}=& \underset{\mathbf{x}\in X(t_{\max })}{\max }{f}_{\text{cost}}(\mathbf{x})-\underset{\mathbf{x}\in X(t_{\max })}{\min }{f}_{\text{cost}}(\mathbf{x})\text{mbox}and\\ L_{\text{LF}}=& \underset{\mathbf{x}\in X(t_{\max })}{\max }{f}_{\text{LF}}(\mathbf{x})-\underset{\mathbf{x}\in X(t_{\max })}{\min }{f}_{\text{LF}}(\mathbf{x})\end{array}$$

$$\left[\begin{array}[]{cc}\frac{\underset{\bm{x}\in X(t_{\max})}{\min}f_{\text{cost}}(\bm{x})}{L_{\text{cost}}}&\frac{\underset{\bm{x}\in X(t_{\max})}{\max}f_{\text{LF}}(\bm{x})}{L_{\text{LF}}}\\ \end{array}\right]^{T}$$

$$O(2t_{\max }{N}_{\text{nom}}(t_{\text{com}}+N_{\text{nom}}-1))$$

Full text

I Introduction

In existing power grids, power plants usually deliver a unidirectional power flow to customers, converting only one-third of the total energy in their fuel into electricity and wasting the heat produced. Twenty percent of a power grid’s generation capacity is often used only to cover peak loads that account for approximately five percent of the time [1, 2]. When natural disasters occur, conventional power grids are unable to resist or self-heal. Next-generation power grids, known as smart grids, have been proposed and designed to replace traditional power grids with the aim of reducing transmission loss, generating electricity more efficiently, and resisting or self-healing after natural disasters. A few governments have adopted proactive policies to popularize smart grids and construct related infrastructure. Smart grids are expected to have higher electricity transmission stability, detect faults and self-heal, allow bidirectional energy and information flow between customers and suppliers, and reliably defend against attacks or natural disasters [3, 4].

In a smart grid, the real-time energy consumption profiles of residential users are crucial to power suppliers[5]. With this information, suppliers or local trading centers can design a dynamic pricing scheme that strengthens market mechanisms and encourages users to shift their peak loads [6, 7]. Power scheduling can then be achieved through the direct minimization of energy consumption costs [8, 9, 10, 11, 12, 13, 14, 15]. The change of demand curves in response to pricing signals is termed demand response.

Common pricing models include real-time pricing, day-ahead pricing, time-of-use pricing, and critical-peak pricing models[16]. Regardless of which dynamic pricing scheme is employed, suppliers can lower the cost of power generation if users’ peak loads are shifted. This can be achieved through the use of demand response programs. Generally, shifting peak loads is related to increasing the load factor, which is the reciprocal of the peak-to-average ratio [17, 18]. However, to fully utilize the grid capacity, the load factor is a more appropriate metric.

Several studies have incorporated the concept of load factor or peak-to-average ratio into power scheduling. Game theory approaches have been widely used for designing demand response programs[19, 20, 21, 22]. Although under certain conditions Nash equilibrium strategies can be Pareto optimal[23], it is not always the case: payoffs of game players may be improved simultaneously after a Nash equilibrium is attained[24]. Stochastic optimization such as genetic algorithms has been applied to attain a satisfactory load factor as well[25], yielding single-objective optimization problems. Some drawbacks may inherit from such single-objective formulations, such as the need of prior decision making for the tradeoff between objectives[26].

Because demand response is closely related to pricing signals, a demand response program can be derived from a pricing scheme[27]. Kunwar et al.[28] proposed an area-load based pricing scheme for demand side management. The load factor and energy cost were jointly optimized. Given pricing signals, the proposed methodology could produce a promising demand response program addressing the two objectives. In contrast with single-objective optimization, such a multiobjective approach could avoid, for example, the need of heuristic assignment for weighting coefficients and the need of prior decision making for tradeoffs induced by conflicting objectives. There are, however, a few important issues that were not covered by[28]. First, although the approach considered shiftable loads, they were addressed statistically. In this case, no integer or discrete decision variables were involved in optimization, but they are important for explicit control of shiftable loads. Second, owing to the employed statistical modelling, the impact of electric vehicles (EVs) was not explicitly examined. In general, EVs pose different physical constraints and should be distinguished from ordinary home appliances. Third, renewable energy sources (RESs) that have been widely adopted in residential communities were not considered.

Motivated by[28], we propose a multiobjective approach that addresses the load factor and energy consumption costs of residential users in separate dimensions, leading to a multiobjective optimization problem (MOP). Our system models include EVs and RESs. Regarding the MOP, power scheduling profiles serve as the decision variables, consisting of continuous and discrete ones. For these two types of decision variables, feasible search mechanisms are developed by exploiting constraint and variable structures. Relevant analysis is provided to show their effectiveness. Solving the MOP yields an approximate Pareto set and a Pareto front. A Pareto optimal solution is then selected to indicate how power loads should be adjusted over time, yielding a Pareto optimal demand response (PODR) program. This program can be offered by utilities to residential users. The associated algorithms or technologies are implemented in residential homes. The load factor is related to grid reliability for utilities while the energy cost pertains to the participation of residential users.

The main contributions of this study are summarized as follows: First, we examine the grid load factor and energy costs of residential users using a multiobjective framework, and consider a combination of various types of appliances, an EV, and RESs. This combination has not been fully investigated in the literature when multiobjective optimization is applied for residential power scheduling. Second, we propose a few algorithms that generate feasible values for decision variables, leading to a solution method for our MOP. Owing to physical constraints involving discrete and continuous variables, this MOP can hardly be solved by conventional multiobjective evolutionary algorithms. Third, numerical simulations demonstrate that our approach can be beneficial to residential users while achieving a satisfying load factor as compared to comparable methods. Finally, in addition to numerical evaluations, the proposed algorithms are theoretically analyzed, providing a rigorous foundation for the proposed PODR program.

The remainder of this paper is organized as follows: Section II discusses related work. Section III describes mathematical models of various home appliances, the energy consumption costs of residential users, and the load factor of the power grid. The PODR program is proposed in Section III. Section IV presents our simulation results involving comparisons between existing demand response strategies. Finally, the paper concludes in Section V.

II Related Work

This section examines various methods employed to design demand response programs for residential users. To facilitate discussions, the section is divided into two subsections. The first subsection presents existing residential demand response methods involving the concept of the load factor but without explicitly addressing it as an objective. This includes methods that have the load factor as an optimization constraint or have the load factor numerically analyzed in the simulations. The second subsection investigates methods that do not involve the concept of the load factor but pertain to our theme.

II-A Methods With Load Factor

Game theory, control methods, mixed integer nonlinear programming, convex optimization, and machine learning methods have been the most promising techniques employed to realize residential demand response involving the concept of the load factor. When game theory is used to model demand response programs[29, 30, 31, 32, 33, 34, 35], two games are typically involved: one for suppliers such as utilities and aggregators, and the other for customers. A demand response program is realized by attaining the Nash equilibrium of the games. Game theory approaches often involve bidding, but a bidding scheme for residential demand response can also be realized using other techniques, see, for example, [36]. Game theory approaches can be further combined with blockchain technologies to improve network security[37].

Control strategies have been applied for demand response services as well[38]. In this case, state-space representations are mostly employed to model system dynamics. For instance, Luo et al.[39] proposed a multistage home energy management system consisting of forecasting, day-ahead scheduling, and actual operation. At the day-ahead scheduling stage, a coordinated home energy resource scheduling model constrained by the peak-to-average ratio was constructed to minimize a one-day home operation cost. At the actual operation stage, a model predictive control based operational strategy was proposed to correct home energy resource operations with the update of real-time information.

Models for residential energy demand have drawn much attention[40, 41]. When models are certain or can be reliably estimated, mixed integer nonlinear programming that involves discrete decision variables can be applied[42, 43, 44, 45]. In this case, discrete or integer decision variables often represent shiftable loads that can be adjusted to change demand curves. If associated problem formulations do not involve discrete or integer decision variables, then conventional linear programming is applicable to residential demand management[46].

Convex optimization methods have been applied to solve a subproblem for demand response programs[47, 48]. In[47], for example, an optimization problem that addressed a trade-off between payments and the discomfort of users was formulated, yielding both integer and continuous variables. The problem was further separated into two subproblems, and one of them was solved using convex optimization methods.

To relax assumptions on entities in the system or probably to address system uncertainties, reinforcement learning for demand response has been extensively studied recently[49]. Through the learning process, an agent or agents can learn how to optimize users’ consumption patterns. Considering the peak-to-average ratio, Bahrami* et al.[50] proposed an actor-critic structure that reduced the expected cost of users in the aggregate load. Dehghanpour et al.*[51] considered air conditioning loads as agents and optimized their consumption patterns through modifying the temperature set-points of the devices. Both consumption costs and users’ comfort preferences were addressed.

Given the aforementioned state-of-the-art methodologies, a few points should be noted. While shifting loads through demand response programs proves promising, it is worth mentioning that load shifting may have impact on distribution systems[52, 53]. In addition to routine operations, demand response programs can be utilized for emergency operations[54]. For a comprehensive review, the reader can refer to[55] and[56] on residential demand and to[57] on a case study illustrating the benefits of demand response in consideration of residential users and utilities’ load factor.

II-B Methods Without Load Factor

This subsection examines recent studies on demand response programs without explicitly involving the concept of the load factor; an energy management system has been a dominant way of realizing a demand response program in response to utility pricing signals in the literature. Here is a list of examples. Hansen et al.[58] investigated a home energy management system that automated the energy usage. Observable Markov decision process approaches were proposed to minimize the household electricity bill. Shafie-Khah et al.[59] proposed a stochastic energy management system that considered uncertainties of the distributed renewable resources and the EV availability for charging and discharging. Adika and Wang[36] examined a day-ahead demand-side bidding approach that maximized residential users’ benefits. Rastegar et al.[60] constructed a two-level framework for residential energy management. Customers minimized their payment costs and sent out the desired power scheduling of appliances at the first level; a multiobjective optimization framework was employed to improve technical characteristics of the distribution system at the second level.

Some authors focuses on the reduction of peak loads in residential demand response. Vivekananthan et al.[61] investigated a reward based demand response algorithm that could shave network peaks. Zhou et al.[62] established a multiobjective model of time-of-use and stepwise power tariff for residential users, yielding load shifting from peak to off-peak periods.

Most recently, learning based approaches to residential demand response have emerged, see, for example, [63, 64] and [65]. In addition to direct investigations on methodologies, tools have been developed to facilitate existing demand response processes. For instance, Paterakis et al.[66] employed artificial neural networks and the wavelet transform to predict the response of residential loads to price signals. Wang et al.[67] developed a method for estimating the residential demand response baseline.

Finally, it is worth mentioning that residential user behaviors may be studied in an aggregate manner, thereby introducing the concept of residential demand aggregation[68]. Two emerging topics are demand response methods for residential heating, ventilation, and air conditioning[69, 70, 71] and demand response methods for residential plug-in EVs[72, 73, 74].

III System Models

This section discusses mathematical models describing the power consumption of home appliances, energy costs of residential users, and load factor of the grid. The day-ahead pricing scheme is considered. To realize an automatic demand response program, we consider an Internet of Things based environment with home appliances able to transmit and receive signals[75, 76]. These home appliances can be controlled by a central controller using IEEE standards, such as IEEE 802.3, IEEE 802.11, or IEEE 1901. Fig. 1 presents the system diagram and Table I summarizes notations and acronyms used throughout this paper.

III-A Home Appliances

Home appliances can be classified into three types [77]: unshiftable and inflexible appliances (the corresponding set is denoted by $\mathcal{A}$), shiftable and inflexible appliances (the corresponding set is denoted by $\mathcal{A}_{S}$), and shiftable and flexible appliances (the corresponding set is denoted by $\mathcal{A}_{SF}$). The loads induced by using those appliances can be classified accordingly. Let $h\in{\mathcal{H}}$ denote the time slot and $\Delta h$ denote the slot duration.

III-A1 Unshiftable and Inflexible Appliances

Some home appliances are used during specific time periods. For example, lights must be turned on in the evening, and it may not be possible to shift the time of use or adjust their switching time. Cooking appliances, such as electric pots, roasters, and microwave ovens, are also used in specific time slots. Residential users may use entertainment electronics such as televisions and computers only after work or school. These appliances are thus regarded as unshiftable and inflexible. For $a\in{\mathcal{A}}$, we denote $p_{a}^{h}$ as the amount of power for appliance $a$ to operate in time slot $h$.

III-A2 Shiftable and Inflexible Appliances

Shiftable and inflexible appliances are used at potentially any time, but their power consumption when performing a specified job is fixed (and thus inflexible). Vacuum cleaners and washing machines belong to this type. The time at which an automatic vacuum machine or a washing machine is operated may not be crucial (assuming that they do not make too much noise; otherwise, they should not be operated at night). A user could prescribe a possible time window, and then the demand response program could decide when the device operates according to these instructions. The time window must be sufficient for the appliance to complete the job.

Let $w_{b}$ denote the total working hours required by appliance $b\in{\mathcal{A}_{S}}$. Suppose that a residential user sets the acceptable start and end time slots of appliance $b$’s operation as $s_{b}$ and $e_{b}$, respectively ($s_{b}\leq e_{b}$). Let ${C_{w_{b}}^{s_{b}\to e_{b}}}$ denote the set of all $w_{b}$-combinations of working time slots within the time window $[s_{b},e_{b}]$. A demand response program thus selects certain working time slots represented by $h_{b}$ from ${C_{w_{b}}^{s_{b}\to e_{b}}}$ for appliance $b$, i.e.,

[TABLE]

For instance, if $s_{b}=15$, $e_{b}=17$ , and $w_{b}$ = 2 (i.e., appliance $b$ needs 2 hours to finish its work), then

[TABLE]

and $h_{b}$ could be equal to $\left\{15,16\right\}$, $\left\{15,17\right\}$ or $\left\{16,17\right\}$. If $h_{b}=\left\{15,16\right\}$, then appliance $b$ operates from 14:00 to 16:00. The associated power for appliance $b$ to operate in time slot $h$ is denoted

[TABLE]

where $p_{b}^{h}(h_{b})=0$ if $h_{b}\cap\left\{h\right\}=\emptyset$.

III-A3 Shiftable and Flexible Appliances

The use of shiftable and flexible appliances can be adjusted in two dimensions: when the appliances are used and how much power they should consume. Heaters, air conditioners, and EV charging stations can yield shiftable and flexible loads. For example, residential users can lower an air conditioner’s temperature setting (flexibility) and delay its use (shiftability) as they wish [9].

Let $p_{c}^{h}$ denote the amount of power for appliance $c$ to operate in time slot $h$, where $c\in\mathcal{A}_{SF}$, and $s_{c}$ and $e_{c}$ denote the start and end time slots of appliance $c$, respectively. The following constraints are imposed on $p_{c}^{h}$ [13]:

[TABLE]

where $p_{c}^{\min}$ and $p_{c}^{\max}$ represent the minimum and maximum operating power, respectively, and $p_{c,\text{agr}}^{\min}$ represents the aggregate minimum power for using appliance $c$, which pertains to user comfort. To ensure that feasible $p_{c}^{h},h=s_{c},...,e_{c},$ exist, we assume

[TABLE]

Although charging an EV can be considered as a shiftable and flexible load, we excluded EVs from $\mathcal{A}_{SF}$ because the charing power of an EV depends on the remaining capacity of the vehicle’s battery, which imposes additional constraints [78]. Let $p_{\text{ev}}^{h}$ denote the EV charging power in time slot $h$, and $s_{\text{ev}}$ and $e_{\text{ev}}$ be the start and end time slots of the battery charging, respectively. The following constraints on the charging rate and EV battery capacity must be satisfied:

[TABLE]

where $p_{\text{ev}}^{\max}$ represents the maximum charging rate, $B_{\text{ev}}^{\min}$ is the minimum capacity of the EV battery, $B_{\text{ev}}^{0}$ is the initial capacity, and $B_{\text{ev}}^{\max}$ is the maximum capacity. To ensure that feasible $p_{\text{ev}}^{h},h=s_{\text{ev}},...,e_{\text{ev}},$ exist, we assume

[TABLE]

III-B Energy Cost and Load Factor

A residential home can be integrated with renewable energy sources (RESs) and equipped with a storage system for energy management. Let $B^{h}$ be the energy level of the storage system ($B^{0}$ represents the initial energy level), $r^{h}$ be the expected charging power from RESs, and $u^{h}$ be the control law that dictates the amount of power being charged to or discharged from the energy storage system. The storage dynamics can be described as[79]

[TABLE]

The constraints on the energy storage system are

[TABLE]

where $B^{\max}$ is the maximum storage capacity.

The total energy extracted from the grid in time slot $h$, denoted by $E_{\text{total}}^{h}$, can be expressed as

[TABLE]

If $u^{h}$ is greater than the total power demand, then $E_{\text{total}}^{h}=0$ and the excess energy is discarded. This situation can happen when $r^{h}$ is too large to be stored in the energy storage system and consumed by residential appliances. If $u^{h}<0$, then $-u^{h}$ is the amount of power delivered to the energy storage system from the power grid. The total energy cost can be obtained by multiplying the total energy consumption by the electricity price $\lambda^{h}$:

[TABLE]

Unlike customers, utilities are principally concerned with the load factor. The load factor is critical because generation cost and grid quality have a direct connection with the load factor. A higher load factor implies a more stable power grid and a lower cost of power generation[80]. The load factor can be defined as the ratio of the average demand to the maximum demand [17, 18]. The following performance index can be used in a demand response program offered by utilities to improve their load factor:

[TABLE]

where $H$ is the size of $\mathcal{H}$ and represents the observation horizon.

IV Pareto Optimal Demand Response Program

This section proposes the PODR program that can be offered by utilities to residential users. An MOP pertaining to power scheduling is formulated. The objectives of the optimization problem are to minimize the energy cost of a residential user in (9) and maximize the load factor in (10). To construct a stochastic search scheme for exploration and exploitation of the decision space, we analyze the associated decision variables and propose a few methods of feasible value generations and mutation and crossover operations that preserve feasibility. A multiobjective evolutionary algorithm for constructing the PODR program is developed accordingly. Finally, the algorithm complexity is discussed.

The MOP is formulated as

[TABLE]

where $\bm{x}$ denotes the decision variable vector that contains decision variables $h_{b}$, $p_{c}^{h}$, $p_{\text{ev}}^{h}$, and $u^{h}$ for all $b,c,$ and $h$. The objective functions are conflicting, so the global optimal solution that optimizes both objective functions simultaneously does not exist. To solve the MOP in (11), Pareto optimality is adopted[81]. A feasible point $\bm{x}^{\prime}$, i.e., an $\bm{x}^{\prime}$ that satisfies the constraints, dominates another feasible point $\bm{x}^{\prime\prime}$ if the conditions $f_{\text{cost}}(\bm{x}^{\prime})\leq f_{\text{cost}}(\bm{x}^{\prime\prime})$ and $f_{\text{LF}}(\bm{x}^{\prime})\geq f_{\text{LF}}(\bm{x}^{\prime\prime})$ hold true with at least one strict inequality. A solution is nondominated (also called Pareto optimal or Pareto efficient) if improving one objective value must yield a degradation in the other objective value[82]. A set of Pareto optimal solutions or nondominated solutions is desired. A nondominated solution should be selected from the set on the basis of its associated performance represented by the Pareto front.

In the following discussions, we adopt a few terminologies used in genetic algorithms[83]. Genetic algorithms are stochastic search techniques that have roots in genetics and employ a population based method to find solutions to optimization problems conventionally with a single objective. These algorithms begin with a set of points randomly initialized. The set is termed the initial population. Points in the population are evaluated on the basis of the objective function, yielding function values called the fitness. With the help of the fitness, a set of new points are generated using mutation and crossover operations. These operations together with a selection operation based on the fitness are applied to improve an average fitness value from population to population. The aforementioned procedure repeats iteratively to produce new populations until a stopping criterion is satisfied.

In genetic algorithms, the mutation is performed on a candidate point and derives a new point termed a mutant. The probability of having a mutant is dictated by a mutation rate. If the candidate point is represented by a binary string using an encoding scheme, then a typical mutation can be designed as stochastically complementing bits from 0 to 1 or vice versa. While the mutation is applied to one candidate point, the crossover is performed on a pair of candidate points called the parents and produces a corresponding pair of points called the offspring. The probability of performing the crossover operation is dictated by a crossover rate. In the case of using the binary representation, a typical crossover operation can be realized through an exchange of substrings of the parents.

To address constraints of optimization problems in genetic algorithms, penalty methods can be used. A penalty function is added to the objective function for point evaluation. The fitness of a point becomes a sum of the objective function value and a penalty function value. For an infeasible point, i.e., violating the constraints, its penalty function value is nonzero (negative for maximization problems and positive for minimization problems) and thus penalizes the objective value. Because points with better fitness are prone to be kept in populations during the algorithm iterations, infeasible points can be gradually removed, leading to a feasible set.

Genetic algorithms provide a generic framework for solving optimization problems. Suitable modifications can render them more efficient and powerful. For example, we may use certain schemes dedicated to the problem of interest that randomly generate feasible points and provide mutation and crossover operations that preserve the feasibility of those points. As such, algorithm efficiency can be improved because more computations are spent on improving feasible points instead of finding feasible ones. Furthermore, by incorporating the concept of Pareto optimality into the fitness evaluation, we may design algorithms that can address multiple objectives, which are termed multiobjective evolutionary algorithms.

Although a few multiobjective evolutionary algorithms are available for finding solutions of MOPs, most of them consider either pure continuous decision variables or pure discrete decision variables. The situation in which both continuous and discrete decision variables are involved, as in our scenario, has not been well addressed[84]. In addition, the constraints in (11) can be difficult to address using conventional constraint handling techniques. Methods of feasible value generations and mutation and crossover operations that preserve feasibility are required.

To solve (11), we first consider an encoding scheme for discrete variables and present the associated mutation and crossover operations. Continuous variables are then addressed. For discrete variable $h_{b}$, a feasible value can be readily generated according to (1). The following encoding scheme and associated mutation and crossover operations are adopted. Let $\phi(h_{b})$ be a binary representation of $h_{b}$ with length $e_{b}-s_{b}+1$ and $[\phi(h_{b})]_{j}$ denote the $j$th bit of $\phi(h_{b})$. The encoding scheme $\phi(\cdot)$ is designed as

[TABLE]

for $j=1,2,...,e_{b}-s_{b}+1$. For instance, if $s_{b}=15,e_{b}=17,w_{b}=2$, and $h_{b}=\left\{15,16\right\}$, then $\phi(h_{b})=\;$110; if $h_{b}=\left\{15,17\right\}$, then $\phi(h_{b})=\;$101. For the mutation operation, we randomly find a pair $([\phi(h_{b})]_{j},[\phi(h_{b})]_{j^{\prime}})=(0,1)$ and switch their values. For the crossover operation, we apply logic operation “OR” to two binary representations, and then randomly choose $w_{b}$ bits with value 1 from the offspring and convert the other bits with value 1 to value 0. The condition in (1) is thus satisfied through the use of the mutation and crossover operations.

For continuous decision variables $p_{c}^{h},p_{\text{ev}}^{h},$ and $u^{h}$, we note that the following mutation and crossover operations can produce feasible points in the decision space if points involved are all feasible.

Theorem 1

Let $x_{\text{new}}^{h}$ denote a mutant or offspring. The mutation or crossover is performed according to

[TABLE]

where $\delta\in[0,1]$ is a random number for all $h$. For the mutation operation, $x^{h}$ is a point in the population and ${x^{\prime}}^{h}$ is a point that is randomly generated. For the crossover operation, $x^{h}$ and ${x^{\prime}}^{h}$ are distinct points in the population. If $x^{h}$ and ${x^{\prime}}^{h}$ in (13) are feasible, then the mutant or offspring $x_{\text{new}}^{h}$ is feasible.

**Proof: **Consider $x^{h}=p_{c}^{h}$. If $p_{c}^{h}$ and ${p^{\prime}}_{c}^{h}$ are feasible, then

[TABLE]

We have

[TABLE]

for $\delta\in[0,1]$. Therefore, $p_{\text{new}}^{h}=\delta p_{c}^{h}+(1-\delta){p^{\prime}}_{c}^{h}$ is feasible. A similar argument can be performed to show the feasibility property when $x^{h}=p_{\text{ev}}^{h}$.

Consider $x^{h}=u^{h}$. If $u^{h}$ and ${u^{\prime}}^{h}$ are feasible, then we have

[TABLE]

Define

[TABLE]

By (14), we have

[TABLE]

Note that

[TABLE]

Therefore, $B_{\text{new}}^{h}\in[0,B^{\max}]$, which implies that $u_{\text{new}}^{h}$ is feasible. $\Box$

If there are methods of randomly generating feasible points, then the feasibility can be preserved during the solution process using the mutation and crossover operations in Theorem 1. Algorithms 1, 2, and 3 present those methods, and their effectiveness are further confirmed by the following theorem.

Theorem 2

Under the assumptions described in (3) and (5), Algorithms 1, 2, and 3 produce feasible values for $p_{c}^{h},p_{\text{ev}}^{h},$ and $u^{h}$.

**Proof: **According to (15) and the fact that $p_{c}^{h},h=s_{c},...,e_{c},$ are initially generated from $[p_{c}^{\min},p_{c}^{\max}]$, variables $p_{c}^{h},h=s_{c},...,e_{c}$, with $p_{c}^{h}\geq p_{c}^{\min}$ approach $p_{c}^{\max}$ from the left as the number of iterations increases. Owing to the assumption in (3), the conditions in (2) hold true eventually when the values of $p_{c}^{h},h=s_{c},...,e_{c},$ increase. By a similar argument and according to (5) and (16), we have $B_{\text{ev}}^{\min}\leq B_{\text{ev}}^{0}+\sum_{h=s_{\text{ev}}}^{e_{\text{ev}}}p_{\text{ev}}^{h}\Delta h$ eventually when the values of $p_{\text{ev}}^{h},h=s_{\text{ev}},...,e_{\text{ev}},$ increase; however, it is possible that an increment is too large to have $B_{\text{ev}}^{0}+\sum_{h=s_{\text{ev}}}^{e_{\text{ev}}}p_{\text{ev}}^{h}\Delta h\leq B_{\text{ev}}^{\max}$. To remedy this, we use (17) to reduce the values of $p_{\text{ev}}^{h},h=s_{\text{ev}},...,e_{\text{ev}}$. When (17) is executed, we have the following two results: $p_{\text{ev}}^{h}$ on the left-hand side of (17) satisfies $p_{\text{ev}}^{h}<p_{\text{ev}}^{\max}$ because the term $(B_{\text{ev}}^{\max}-B_{\text{ev}}^{0})/\sum_{h=s_{\text{ev}}}^{e_{\text{ev}}}p_{\text{ev}}^{h}\Delta h$ on the right-hand side of (17) is less than 1; and $p_{\text{ev}}^{h}$ on the left-hand side of (17) satisfies $\sum_{h=s_{\text{ev}}}^{e_{\text{ev}}}p_{\text{ev}}^{h}=(B_{\text{ev}}^{\max}-B_{\text{ev}}^{0})/\Delta h.$ The conditions in (4) are then satisfied. Finally, we note that (18) implies (7), illustrating the feasibility of $u^{h}$. $\Box$

With the help of Algorithms 1–3, Algorithm 4 presents the PODR program that can be offered by utilities to residential users.

In Algorithm 4, information about appliances and acceptable times of use is set first. Mutation and crossover operations are then performed over $X(t_{c})$ in Step 2.1. In Step 2.2, the archive is updated by removing some nondominated points if the population size $\mid X(t_{c})\mid$ is too large. After a number of iterations, an approximate Pareto set and front are obtained in Step 3. In practice, the knee of the approximate Pareto front is often preferred for several reasons: it can achieve excellent overall system performance if the front is bent; it represents the solution closest to the ideal one that is not reachable; and it has rich geometrical and physical meanings[85, 86]. In Step 3, the knee solution is selected according to [87]:

[TABLE]

where

[TABLE]

are the maximum spreads of the approximate Pareto front in the first and second dimensions, respectively, and

[TABLE]

is the ideal vector.

As shown in (19), the final Pareto optimal solution is selected using the system’s knowledge of the maximum spreads, critical information in multicriteria decision making. Such information is not obtained using conventional single-objective optimization approaches. The knee solution $\bm{x}^{*}$ corresponds to a Pareto optimal power consumption profile, yielding the PODR program.

The complexity of Algorithm 4 mainly depends on the use of Pareto concepts that are the most computationally expensive. The complexity can be roughly expressed as[84]

[TABLE]

where $t_{\text{com}}$ represents the average computational time for objective evaluation.

V Simulation Results

This section presents an examination of the power scheduling of 400 residential users.111In [18], the total available capacity of a region was set at 2296 kW. In our scenario, each residential user consumes 5.5 kW at most, which thus accounts for approximately 400 residential users in the region. Each residence was equipped with at most thirteen unshiftable and inflexible appliances, three shiftable and inflexible appliances, two shiftable and flexible appliances, one EV, and possibly one energy storage system integrated with solar energy sources. The exact number of appliances and associated types for a residential user were randomly chosen. Solar energy data in[88] were used. Let ${\mathcal{H}={\left\{1,2,...,24\right\}}}$ and $\Delta h=1$. Table II presents the associated settings, and the notation “unif($h_{1},h_{2},t$)” therein indicates that an appliance required $t$ hours of working time starting from a number randomly chosen from $\{h_{1},h_{1}+1,...,h_{2}\}$. For example, appliance $a_{4}$ needed 5 hours to complete the task, and the start time slot could have been 16, 17, or 18. Most parameter values were chosen according to[89] and[90].

The day-ahead prices illustrated in Fig. 2 were used [91]. After analyzing total energy consumption of each month in 2017, we discovered that energy consumption peaked in July. Four representative days in July were selected and listed in a decreasing order in terms of their price range as follows: 27 July 2017 (maximum price range), 6 July 2017 (25th percentile), 23 July 2017 (75th percentile), and 29 July 2017 (minimum price range). Among them, 27 July 2017 also had the highest price.

The PODR program for users was obtained using Algorithm 4 with $\mu=0.8,N_{\text{nom}}=40,N_{\max}=400$ and $t_{\max}=400$. Approximate Pareto fronts associated with residential users were produced. Fig. 3 plots a sampled approximate Pareto front. Each point on the front corresponded to two values: energy cost and load factor. The shape of the front confirmed that minimizing the cost and maximizing the load factor were conflicting objectives. The knee was selected according to (19).

Our multiobjective approach was compared with an area-load method modified from[28], a payment minimization method modified from[18], and load factor maximization and load variance minimization methods modified from[92] (see Appendix). Owing to constraint structures and different types of variables that were involved, conventional stochastic search methods could hardly produce feasible points. Algorithms 1–3 were thus applied to produce a portion of initial points that were used through the solution processes associated with the comparable methods. Fig. 4 presents individual performance on the representative days. Table III summarizes the performance (expressed as a percentage) of each method. For the energy cost, the percentages were evaluated with respect to our method. A smaller percentage indicated a lower cost. For the load factor, the percentages were evaluated with respect to the load variance minimization method. A larger percentage indicated a higher load factor.

Among those methods, the PODR program balanced the two conflicting objectives in an advantageous manner. This should be generally the case because our approach considered join optimization of the two objectives, had robust constraint handling techniques presented in Algorithms 1-3, and employed dedicate mutation and crossover operations to preserve solution feasibility, as shown in Theorem 2. The modified area-load method also adopted a multiobjective approach; as compared to the PODR, it yielded a small improvement in the load factor by approximately 2.8% (-7.8%+10.6%) but a large degradation in the energy costs by 14.5%. For single-objective optimization methods, load variance minimization and load factor maximization yielded larger load factors than other methods because performance metrics related to the load factor were optimized; however, the resulting energy costs for residential users, which were not considered during the solution process, were the worst among other methods. Finally, the modified payment minimization method yielded satisfactory costs for residential users but the worst load factor; the associated energy cost was higher than that of the PODR because that method lacked suitable constraint handling techniques.

VI Conclusion

This paper investigated two critical aspects of power scheduling: residential users’ energy costs and the power grid’s load factor. These two aspects should be considered simultaneously when utilities are to offer demand response programs to residential users. In practice, minimizing the cost and maximizing the load factor are two conflicting objectives—one cannot be optimized without worsening the other. To address this challenge, we proposed a Pareto optimal demand response program that determined the Pareto optimal power scheduling. This program was constructed by solving a multiobjective optimization problem. To explore and exploit the decision space associated with the problem, a few algorithms were developed to generate feasible values for decision variables. Relevant analysis was performed to show the effectiveness of the algorithms. Compared with existing methods for power scheduling, the Pareto optimal demand response program found a satisfying balance between the two objectives by sacrificing one objective slightly while substantially improving the other. This study was mainly focused on demand response for residential users. Our future work includes the investigation on multiobjective approaches to the design of demand response programs for commercial or industrial customers and for vehicle-to-grid systems.

[Descriptions of Comparable Methods]

For modified area-load and payment minimization methods, the constraints in (11) were addressed using the following constraint function:

[TABLE]

If $U(\bm{x})>0$, then the power scheduling profile $\bm{x}$ violates the constraints and $\bm{x}$ is infeasible.

For the area-load method modified from[28], we solved

[TABLE]

where $\sigma_{1}$ and $\sigma_{2}$ represent the weighting factors of the constraint violation and penalty, respectively. Constraint handling was not discussed in [28], but it is essential because physical systems always induce constraints. We modified the objectives by including $U(\bm{x})$. The penalty of energy load deviating from the average load was evaluated as

[TABLE]

where

[TABLE]

The nondominated sorting genetic algorithm-II was used to solve (21).

For the payment minimization method modified from[18], we solved

[TABLE]

using particle swarm optimization.

The load factor maximization method from[92] was modified as

[TABLE]

For the load variance minimization method modified from[92], we solved

[TABLE]

where $E_{\text{av}}$ is defined in (23).