A Closed-loop Control Strategy for Air Conditioning Loads to Participate in Demand Response

Thermostatically controlled loads (TCLs), such as air conditioners (ACs), are important demand response resources—they have a certain heat storage capacity. A change in the operating status of an air conditioner in a small range will not noticeably affect the users' comfort level. Load control of TCLs is considered to be equivalent to a power plant of the same capacity in effect, and it can significantly reduce the system pressure to peak load shift. The thermodynamic model of air conditioning can be used to study the aggregate power of a number of ACs that respond to the step signal of a temperature set point. This paper analyzes the influence of the parameters of each AC in the group to the indoor temperature and the total load, and derives a simplified control model based on the two order linear time invariant transfer function. Then, the stability of the model and designs its Proportional-Integral-Differential (PID) controller based on the particle swarm optimization (PSO) algorithm is also studied. The case study presented in this paper simulates both scenarios of constant ambient temperature and changing ambient temperature to verify the proposed transfer function model and control strategy can closely track the reference peak load shifting curves. The study also demonstrates minimal changes in the indoor temperature and the users' comfort level.


Introduction
Demand Response (DR) is defined as the changes in electric usage by customers from their normal usage patterns with changes in the price of electricity over time.It is considered a critical application in smart grids.With the Advanced Metering Infrastructure (AMI), the power usage of different appliances can be adjusted either directly, e.g., through incentive-based programs (IBP), such as operational parameters/states changing requested by grid operators, direct load control, and interruptible loads; or indirectly, e.g., through price-based programs (PBP) [1], such as real-time pricing, time of use pricing, and critical peak pricing [2,3].By smoothing out the system power demand over time, DR is capable of providing peak shaving, load shifting and ancillary services to achieve system reliability and stability.The performance of DR programs is measured by peak load reduction and demand elasticity.
Traditionally, the power system adopts the strategy of "electricity production determined by loads" to achieve power balance.Loads have been regarded as the passive physics terminal.In recent years, the smart grid concept has become popular in the electric system.The communication capability of the new intelligent service network and the loads' controllability have been greatly improved.Direct Load Control (DLC), which usually reduces the loads by controlling the thermostatically controlled loads (TCLs) of air conditioners, fridges, water heaters of the residential users or the small-business users, is an important type of incentive-based demand response.These kinds of loads are called flexible loads.They have heat energy storage ability and can switch or change the control parameters in a short time, which will not impact users.Meanwhile, they have the potential to balance system power.Flexible loads can have a real-time response to the grid's demand, reduce the backup capacity properly, and improve the level of the power system to run safely and economically [4].They can be incorporated into the normal operation of power system's dispatch by the demand response [5].
Electricity shortages usually occur during peak hours in the summer or winter, when the percentage of the air-conditioning loads reach 30%~40% of the total capacity-the demand peak rises with the number of air conditioners in use, while during off-peak hours, there is usually a power surplus.The types of the air conditioners include centralized and non-centralized [6,7].A centralized air conditioner controls the temperature in different rooms by a mainframe connecting to several terminals through the ventilation system, while non-centralized air conditioners (ACs), which account for a large proportion of all ACs, are installed in individual rooms and control the temperature independently.TCL presents a huge opportunity in the power grid.The effect of controlling the loads of TCLs is equivalent to building a power plant with the same capacity.However, the investment of the demand response is much lower compared with the construction of a new power plant, and reasonable load control has minimal impact on users.The cost for building a peaking power plant is 1200 dollars per kilowatt [8], while the compensation for users about load management is much lower and can be put into use in a short time.
There are several ways to model TCLs, including modeling based on the actual physical process, regression modeling based on historical data, and modeling based on the Fokker-Planck diffusion equations.Other authors focus on a probability characteristics model and the use of black-box model identification techniques [9][10][11].These models are often too complicated for the control system.For example, the physical process model, which needs to be converted to the mathematical model of control, cannot be directly used in numerical simulations.The regression model requires a large amount of historical data.Similarly, the model based on Fokker-Planck equation is very complex and difficult to implement [12].A state sequence model, proposed by Lu and Chassin, considers the uncertainty modeling of thermostatically controlled loads.The disadvantage of this model is that it can only slowly change the load over time.
In most cases, the control scheme of TCLs is based on the idea of optimization.A scheduling decision model has an objective function to maximize the load aggregate interest and to minimize the actual output deviation of ACs.It includes a constraint condition that the indoor temperature has to be within the range of [ min max , T T ] in order to ensure the user's comfort level [13].After the air conditioning model is determined, the air conditioning scheduling scheme is calculated using the optimization algorithm.This method belongs to the open-loop control, with intensive calculation, cannot achieve accurate control and real-time control which is based on load output change.The research by Na et al., on the direct load control (DLC) and controlling AC loads by electricity price in order to achieve the purpose of load reduction does not include control strategies.Manichaikul and Schweppe used the predictive method to control centralized TCLs and provide the auxiliary service for the power grid, taking the minimum error as the target [14].Centralized TCL load management is carried out by using the comfort control method proposed in Malhamé's and Chong's paper.The reference signal which is compared with the actual signal, is the desired aggregated air-conditioning load.Then, the temperature set points were calculated using the intima controller [15].In [16], the electric heat load DLC algorithm is based on the state sequence control algorithm, which aims to stabilize tie-line power fluctuations between distribution network and micro-grid including clean energy.
These methods adopt price signals or power surge signals to control the switch state of the users' equipment.Although they can take part in the ancillary service of the power system, their time scale is usually counted in hours or days, so they cannot follow the target load curve closely.Currently, most load control strategies destroy the diversity of the load group, resulting in power demand oscillation, and produce negative effects to the power system [17].Demand side customers' satisfaction can be taken into consideration using the load management (LM).A fuzzy membership function is established to characterize the users' satisfaction according to the user's continuous control time [18], but it is difficult to express the users' comfort level directly.In this paper, a closed-loop control strategy to solve the above-mentioned problems in existing methods is proposed.
The core contribution and innovation of this paper lies in exploring a simplified two order linear time invariant transfer function model by analyzing the duty cycle of the air conditioning group, proposing a closed-loop control strategy to provide a new way for flexible loads to participate in demand response, and putting forward two suggestions to solve the problem that parameters of the proposed transfer function model vary with the ambient temperature in the case study.Based on the proposed model, accurately load control could be achieved without negatively affecting user experience.
The rest of this paper is organized as follows: Section 2 presents an equivalent thermal parameters model of air conditioning.Section 3 analyzes the influence of AC parameters on its output characteristic and explores a simplified transfer function model for ACs by analyzing the duty cycle of the population.Section 4 presents a design of the model's PID controller based on the PSO algorithm.Section 5 demonstrates the simulation results of both scenarios of constant and changing ambient temperature to verify the proposed transfer function model and control strategy.Section 6 summarizes the discussions and puts forth our conclusions.

An Equivalent Thermal Parameters Model of Air Conditioning
In this section, the correlation between power, temperature and time is established, based on the equivalent thermal parameters (ETP) model.The study focuses on fixed-frequency air conditioning, and uses the Monte Carlo method to simulate the aggregated dynamic response of AC groups.
Air conditioning has a periodic working pattern, as shown in Figure 1: When the air conditioning is on, the temperature inside the house keeps dropping until the temperature reaches the boundary   .Then the air conditioner is off and the temperature inside the house keeps rising until the temperature reaches the boundary   , after which the air conditioning is switched on again.The AC's periodical properties can be described by the following equivalent thermal parameters model [19]: where ( ) is the air temperature inside the i-th house at time t, a  is the outdoor ambient temperature, i C (kWh/°C) is the i-th AC's equivalent thermal capacitance, i R is the i-th AC's equivalent thermal resistance (°C/Kw), i P (kW) is the i-th AC's thermal power, which divided by the coefficient of performance (COP) is the actual load of air conditionings, ( ) i s t is the switch state variable of the i-th AC at time t. ( ) 1 i s t  means the i-th AC is on at time t, while ( ) 0 i s t  means the i-th AC is off at time t.
The formula for calculating the aggregated demand of n ACs in the group is: D t is a per-unit value, taken as the baseline, whose numerator represents the sum of the running AC loads in the group, and the denominator represents the total power of the group, assuming that all ACs in the group are on.
The traditional thermodynamic model of the air conditioning using Equations ( 1)-( 3) has some limitations.Each air conditioning room is a dynamic system with a nonlinear and independent iteration.When there is a large number of ACs, we have to face the problem of dimension disaster which triggers great difficulties.In addition, the model contains both continuous variable temperature θ, and discrete variable switching state s.In the following section, we present a simplified transfer function model, which can be used to easily solve the aggregated load of the population containing n ACs.

Transfer Function Model for Aggregated Air Conditioning Group
It is difficult to perform an accurate mathematical analysis of the traditional thermodynamic model, since each air conditioned room is a dynamic system with a nonlinear and independent iteration.This section explores a simplified transfer function model by analyzing the duty cycle of the air conditioning group.The excitation signal is step change of temperature set point.

The Influence of Air Conditioning Parameters on the Output Characteristic
In order to maintain the diversity of the air conditioning group, we assume that at the initial stage, all air-conditioned room temperatures uniformly distribute between the boundary temperature value [ ,     ], which determines whether the AC is on or off.In order to explore the equivalent heat capacity's relationship with ambient temperature and duty cycle of air conditioning group, we assume that all ACs have the same thermal power P and the same equivalent thermal resistance R.This ensures that the decline rate of the room temperature when the air conditioning is on is equal to the rate of the temperature rise when the air conditioner is off.Therefore, the duty cycle is 0.5.We assume that equivalent thermal capacitance C follows lognormal distribution.Then we have:

2(
) As the control signal, the temperature set point of ACs is the same for all the members in the group: . The hysteresis width that consists of the boundary temperature values is the same for all air conditioners.In the simulation of this paper, the hysteresis width is: Then we analyze the air conditioning group based on the assumption of the above conditions.At t = 0, a step change of 0.5 °C of the temperature set point is occurred.Figure 2a will move right into ( , , which is: Figure 2a,c shows that when the air conditioning is off, the room temperature rises, while Figure 2b,d shows that when the air conditioning is on, the room temperature drops.We define the change rate of the temperature of the i-th air conditioning room as follows: Besides, we define state variables: 0 ( ) where: For an air conditioning group based on the thermodynamic model, at t = 0, a step change of 0.5 °C of the temperature set point is occurred.When ( ) [2 , 2 1]   i x t k k   , the air conditioning is on, and when  , the air conditioning is off.Each time when i x becomes an integer, a on-or-off state transformation occurs.The relationship between i x and the state variable ( ) i s t of air conditioning in the thermodynamic model ( 1)-( 3) is: When ( ) i x t < 1, at time 0+, temperature set point rises.From Figure 3, we can see that only 1/3 of the ACs are on, and by Equation ( 8), the state of the air conditioning is distributed in the interval of ( ) 0.5 i x t  .The probability that an air conditioning in the population switched on is: where Pr[•] is the probability operator.For an air conditioning group based on the ETP model, at t = 0, the temperature set point of all air conditioners rises by 0.5 °C.Assuming that all ACs have the same power, if they are switched on at the same time, we take their total power as the baseline.Therefore, the normalized power demand D(t) of the air conditioning group which consists of n ACs is equal to the probability of a random AC turning on [20]: D(t) also gives the proportion of the air conditioning that is on in the group at time t.The first term of Equation ( 11) is the duty cycle of the air conditioning group at the moment 0+ when the temperature set point rises, while the second term is the duty cycle of the air conditioning group after t > 0. Assuming that all ACs have the same R, P, and the parameter C obeys a lognormal distribution.Consequently, according to Equation ( 6), the temperature change rate v obeys the lognormal distribution, too.We use . .s d to represent the standard deviation, as well as E to the mean value: rel  is defined by: where c  is the mean value of the equivalent heat capacity of n ACs, and c  is the standard deviation.
Based on Equation (11), when obeys the lognormal distribution, it can be considered that ( ) x t is also subject to the same distribution.According to the nature of the lognormal distribution we obtain: where: ( ) ( ) log 1 ( ) From Equation (11) we can also derive: and: The value of the first term of the above two formulas is very small and can be ignored.Therefore: By plugging Equations ( 17) and ( 19) into Equation ( 14), the following equation is obtained: where erf [•] is the Gauss error function.
Let y = 1 and for 2k and 2k + 1, k = 0, 1, … ,.Then, by plugging Equation ( 19) into Equation ( 11), D(t) can be approximated as follows: Figure 4 shows the average temperature and power for the population composed of 10,000 air conditioners according to different values of rel  .At t = 500 min, the temperature set point of all ACs rises 0.5 °C.In Figure 4a-d, for rel  = 0.02, 0.05, 0.2 and 0.5 respectively, the top represents the average temperature inside the house before and after the step response, while the bottom represents the total power of all ACs before and after the step response.rel  means the standard deviation of log-normal distributions as a fraction of the mean value for R, C and P. According to Equation (3), the normalized power demand D(t) means the proportion of the air conditioning which is on in the group at time t.Assuming that all ACs have the same R, P, while C obeys lognormal distribution.This ensures that, the rate of fall of the room temperature when the air conditioning is on, is equal to the rate of rising of the room temperature when the air conditioner is off.The steady-state value of the duty cycle is 0.5.
As shown in Figure 4, the average temperature of the air conditioning group is in the vicinity of the set point 20 °C before a 0.5 °C step.Then, the room temperature stabilizes at 20.5 °C after a period of oscillation.The temperature and power have the same settling time.More generally, the duty cycle may not always remain at 0.5, as in the hotter days when there are more ACs turning on, the duty cycle is larger than 0.5, but in colder days, the situation is the opposite.
When not all parameters such as R, P are the same, for example, as follows a lognormal distribution, variations of temperature and power of the ACs in the group will be more dramatic.Due to the diversity, the damping will be greater.
The temperature and the power of ACs experienced an underdamped oscillation after the temperature set point had risen 0.5 °C at t = 500 min.Figure 4 shows a decaying oscillation.For rel  = / c c   , the bigger rel  is, the greater the system damping becomes.The average temperature and the power curve have a shorter settling time and a smaller overshoot.By observing changes in the aggregated load before and after a step change of the temperature set point, we use a second-order dynamic model to approximate the air conditioning group.It will be easy to solve the aggregated load of a population consisting of n ACs.Thus we can design a controller based on the transfer function later.In the following simulation we can see that, when we release some restrictions on the assumptions in more general cases, the transfer function model and control strategy on the ACs will still be quite effective in real situations.

A Transfer Function Model for the Second Order Linear Time Invariant System
When the temperature set point rises by 0.5 °C, the response of a population consisting of n ACs will be: where ( ) is the aggregated power demand of the AC group given a constant temperature set point.It is also the duty cycle of the steady state.The term including ( ) p G s , which belongs to a transient process of change, indicates the aggregated demand response of a 0.5 °C step in the form of a second-order linear model: where the hysteresis width H is obtained by the upper and lower limit of temperature which determines whether the AC is operating or switched off: The second order transfer function model of an AC group is: Whose parameters are calculated from the equivalent thermal capacitance, thermal resistance, thermal power, hysteresis width and so on.The damping ratio  and the undamped natural oscillation n  of the characteristic polynomials in the denominator are: The coefficients in the numerator are: where v  represents the mean values of v defined in Equation ( 6): where erf represents the Gaussian error function.The proof of the above equation is given in the Appendix A.

Stability Analysis
The aggregated demand response of a 0.5 °C step based on a second-order linear time invariant (LTI) model is shown in Figure 5, which reveals the process of an underdamped oscillation when the system cross the position of equilibrium.
The simulated parameters are listed in Table 1, whose parameters are subsequently calculated by using Equations ( 25)-( 27), according to the air conditioning parameters given in Table 2.  .The decay time of the step response of the control plant to 1% is: ) 520 min.It can be seen that the settling time calculated from the second order LTI model in Figure 5 is consistent with the thermodynamic model in Figure 4. Figure 6 shows the dynamic performance of the open-loop transfer function of the air conditioning group composed of n ACs.We can see from the graph that without the controller, the peak of the two order system is 0.422, the overshoot is 1180%, the rise time is 47.1 s, the adjustment time is 235 s, the system stability is 0.0329, and there exists a static error of 2% for the normalized power demand D(t).Therefore, it is sensible to design a corresponding controller, considering increasing the open loop amplification, and using a proportional control to reduce the static error.

Simulation Analysis
The initial temperature set point of the air conditioning group is assumed to be 20 °C.At t = 500 min, the temperature set point of all ACs is raised by 0.5 °C.The upper and lower limit of the temperature, which determines the operating state of ACs, rises to 20.5 and 19.5 °C.The temperature hysteresis width is still 1 °C.Figure 8 shows a comparison of the aggregated demand between the transfer function model and the equivalent thermal parameters (ETP) model in response to a 0.5 °C step rise.In the ETP model, the number of ACs in the population is 50, 1000 and 10,000 respectively.Figure 8a shows the variation on the average temperature inside the room of the AC group, and Figure 8b shows the variation on the aggregated demand of the AC group.Parameters of the equivalent thermal parameters model can be found in Table 2.When there are only 50 air-conditioning units, although the AC parameters C, R and P are still log-normally distributed in the simulation, the diversity of the population cannot be achieved because of the strong randomness.The simulation curve is very volatile and has many glitches.In contrast, when the number of aggregated ACs in the population is 1000 and 10,000, the population has a strong diversity.The ACs are uniformly distributed in all running states, and the step response curve is very smooth.
We can see that with the parameters in Table 1, the two order linear time invariant transfer function model can be used to approximate the aggregate response of the AC group accurately.The load profile simulated by the transfer function model is very close to which simulated by the ETP model that contains 10,000 ACs.The input signal u(t) represents that, at t = 500, all the ACs' temperature set point has a step rise of 0.5 °C.
The two order LTI model, which has a simple structure, can be used to capture the oscillation characteristics of the aggregate demand response of ACs.This is the key to the design of feedback control strategies for load control.In addition, applicability and complexity can be balanced, which makes the proposed model very practical.
The novelty provided in this section is as follows: the aggregated power of a population of n ACs that respond to the step signal of temperature set point was computed.The dynamics of the response of a numerically simulated population over a realistic range of parameter values were captured.Furthermore, a simplified control model based on the two order linear time invariant transfer function was proposed for the reason that the ETP equations were too complicated to solve.The excitation signal was the step change of temperature set point.The proposed model provided rule-of-thumb about the response with no need for intensive numerical calculations.In addition, it offers facilities for control design.

Proportional-Integral-Differential (PID) Controller Based on Particle Swarm Optimization (PSO) Algorithm
The performance of the PID controller depends on three parameters: Kp, Ki and Kd.Therefore, it is important to optimize these parameters.Currently, the PID controller parameters are mainly adjusted manually, which is not only time consuming, but also cannot guarantee the best performance.The PSO algorithm has been widely used in function optimization, neural network training, pattern classification, fuzzy control system, as well as in other fields [21].This paper will use the PSO algorithm to optimize the parameters of PID controller.Figure 9 shows a load control strategy for the AC group based on the PSO algorithm for PID parameter tuning.The desired reference demand of the ACs minus the actual output power of them is the load tracking error, when the event of load control occurs.This load tracking error is taken as the input of the PID controller, and then, the PSO algorithm is used to adjust its parameters.The calculated temperature set point u(t) is considered as the control signal of the air conditioning group.The bridge between the PSO algorithm and the Simulink model is the corresponding fitness value of the particles, namely, the parameters of the PID controller and the performance index of the system.The optimization process is as follows: firstly, the particle swarm is generated by the PSO algorithm, and the three-dimensional particles in the particle swarm are sequentially assigned to the PID controller parameters Kp, Ki, Kd.Then we run the transfer function model for the AC group, and get the corresponding performance index of the specified parameters.Finally, we take the performance index as the fitness value of PSO to determine whether the exit conditions have been met.If so, the calculation process is terminated, otherwise the particle swarm will continue to update.
The performance index, which is presented by the integral formula of the deviation between the expected system's power demand and the actual feedback output, is a measure of the performance of the control system.In this paper, the instantaneous error function e(t) of the control system is evaluated by the integrated time and absolute error (ITAE) index: The transient response oscillation of the control system based on the ITAE performance index is small.Parameters also have good selectivity.
The novelty provided in this section is as follows: after the analysis of the stability of the proposed model, this paper designs its PID controller based on the PSO algorithm.Compared with the general optimization algorithm, the closed-loop control and feedback strategy are applied to provide a new way for flexible loads like air conditioners to be involved in demand response.
However, the PSO algorithm itself is not the present authors' invention, and it is part of the innovation of this paper.It is applied to tune parameters for the PID controller in order to obtain smaller overshoot, shorter settling time and higher accuracy throughout the control process, and to meet with requires of demand response applications, such as load shifting and peak shaving.The principle and flowchart of PSO algorithm tuning PID controller parameters can be found in Appendix B.

Scenario 1: Simulation under Constant Ambient Temperature Conditions
A population of ACs of similar characteristics can be used as a control group [22,23].We choose the downtown area of the city of Nanjing, China for this research, where there are 18,000 ACs of 10,000 residents participating in the demand response plan.Parameters of ACs in the group are as follows: thermal resistance .When the ambient temperature is constant (e.g., 32 °C), the parameters of the two order transfer function of the aggregate AC group are shown in Table 2.In this scenario the parameters remain unchanged.
First, we do not have any control over the AC population.Therefore the temperature set point would stay fairly stable, and the total power demand of the whole AC population is shown in Figure 10.The temperature of 30 random air-conditioned rooms among the population is presented in Figure 11, and the mean temperature inside the house of all ACs is presented in Figure 12.The simulations were run in MATLAB (The MathWorks, Inc., Natick, MA, USA.) on a desktop computer equipped with an Intel Corei5-3230M 2.60GHz CPU, 4.00 GB memory, and 64 bit Windows 8 operating system.
Figures 10-12 show that, without any control, the AC group will keep running independently.The profile of the aggregated demand of AC group is gentle, the duty cycle in steady state is 0.467.The temperatures of 18,000 air-conditioned rooms are uniformly distributed between 19.5-20.5 °C.The hysteresis width, which determines the upper and lower limits of temperature inside the air-conditioned rooms, is 1 °C.Besides, the mean temperature of all AC rooms is close to the temperature set point 20 °C.
A scenario analysis by varying the temperature set point offset to control the aggregated power demand of air-conditionings is started.Assuming that at the initial stage, 44.6% of the ACs are on, while the rest of them are off, the whole load control scenario is divided into three stages.For the first 100 min, the ACs run independently without any control.The second stage is from 100 to 350 min, during which the grid is in a peak load period.There is a power shortage and therefore needs a 20% reduction of air conditioning load, which is equivalent to 4.3 MW.The load aggregator will raise the temperature set point of the ACs to achieve peak load shifting.The third stage takes place between 350 and 600 min, when the peak load period of the power grid is over, it needs an additional 20% increase of air conditioning load to meet the power consumption of the off-peak period.The load aggregator will drop the temperature set point of the ACs to achieve a satisfactory comfort level.The load management scenario finalizes at t = 600 min.In this scenario, the two order transfer function of the AC group remains unchanged because of the constant ambient temperature, as is shown in Table 2.The particle swarm optimization algorithm is applied to tune parameters for the PID controller.The parameters of the PSO algorithm are set as follows: population size m = 100, dimension D = 3, the 3 parameters Kp, Ki and Kd, to be optimized are in the range of 0-300, the inertia weight w = 0.6, the maximum number of iterations t = 200, acceleration coefficients The IEAT index is chosen as the fitness function, with the minimum fitness value being 0.1.The particle velocity is between [−1, 1].How to select the parameters of the PSO algorithm in detail is explained in [24,25].Figure 13a,b shows the optimization of PID parameters Kp, Ki and Kd by using the PSO algorithm, and Figure 14 shows a variation of the ITAE performance index.With the increase of the number of iterations of the particle swarm, the ITAE performance index, which is considered as the fitness function, gradually stabilizes at around 1.06.Therefore, we get the optimal PID parameters and the performance index as shown in Table 3, in which we tune PID parameters by using the particle swarm optimization algorithm and make a comparison of tuning results with the Ziegler-Nichols method.The parameters ts represents settling time, and tp represents peak time.In order to demonstrate the advantages of the PSO algorithm, we compared it with the classical Ziegler-Nichols method on tuning PID parameters of the two order transfer function model for aggregated ACs [26].Simulation results indicate that, although the peak time of the tuning process is slightly longer, the optimized controller based on the PSO algorithm outperforms the conventional one greatly in the settling time, overshoot and performance index.
Figure 15 shows the comparison between the desired reference value of the aggregated demand of a population of ACs for load shifting, the power output of the population based on conventional PID controller, and the power output of the population based on PSO tuning PID parameters, when the outdoor temperature is a constant 32 °C.In Figure 15, the black curve is the ideal load control reference output, the green curve is the actual power output of the population using the Ziegler-Nichols method to tune PID parameters, and the red curve is the actual power output of the population using the particle swarm optimization algorithm to tune PID parameters.Figure 15b,c is obtained by zooming in on the parts circled by a dotted curve in Figure 15a, and show the effects of reducing the power output in 100-350 min and increasing the power output in 350-600 min, respectively.Table 3. Tuning proportional-integral-differential (PID) parameters using the Ziegler-Nichols method in comparison with the particle swarm optimization algorithm.As can be seen from the figure, using PSO to tune PID parameters is significantly better than using the Z-N algorithm.In response to the step changes of reference values of load control at time 100 and 350 min, the Z-N tuning method has a larger overshoot, a longer settling time, a continuous oscillation throughout the control process and more glitches.In contrast, the control effect on PSO algorithm tuning PID parameters is quite satisfactory, and it meets with requires of rapidity, accuracy and stability of the control system.

Parameters
We can also achieve a 100% load reduction by adjusting the expected reference value of the air conditioning load to zero. Figure 15d shows the simulation results of the aggregated demand of n ACs in this extreme case.We assume that, there appears a serious power shortage and needs a 100% reduction of air conditioning load.The aggregated power consumption of ACs based on PSO tuning PID parameters drops to zero when the reference signal is applied at time t = 100 min, and the whole stage continues for 30 min.
The cost of doing so is that all of the ACs have been switched off and the air temperature inside the house will gradually rise.From a technological perspective, we can achieve any arbitrary adjustment of ACs by using the proposed model and control strategy, but in a real-world scenario, if the duration of 100% load reduction is longer, it will affect the user's comfort and lead to their dissatisfaction on a hot day.Thirty min later, the power demand increases to the initial value, and all of ACs return to the steady state without any parasitic oscillation.
Figure 16 shows the variation of the indoor temperature of 30 random ACs.The temperature hysteresis width H = 1 °C.The air temperature inside the house changes in response to the variation of the temperature set point u(t).We also observe that, the ACs' temperature set point can be controlled by the load aggregator who sets a higher value to reduce the load, as well as setting a lower value to increase the load.After 600 min, the control process is over, but the room temperature does not directly return to the initial value [19.5 °C, 20.5 °C].Instead, it has undergone a nearly 200 min variable process, and the final stable temperature is slightly lower than the initial value.The reason for this phenomenon is that the characteristics of air conditioning itself make the load rebound before and after the load control event.[27][28][29], which reduced two-way broadcasting to one-way and eliminated the problem of unwanted synchronization of TCLs with high quality control and without requiring a channel for consumer-to-grid communication.To quickly compensate for a mismatch between generation and load, the authors put forward strategies based on the timer-based safe protocols (SPs), including SP-1, SP-2 and SP-3, and focus on developing the algorithms to generate a set of power pulses that are useful in spinning reserve applications.
This paper presents a comparison of power demand response between the above hybrid SP protocol combining SP-1 and SP-2 and the proposed closed loop strategy in Figure 17.The simulation scenario is unchanged, from 100 to 350 min, there needs a 20% reduction of air conditioning load (Figure 17a), and an additional 20% increase of air conditioning load to meet the power consumption of the off-peak period (Figure 17b).The green line is the power response under hybrid SP protocol which combines SP-1 (yellow) and SP-2 (blue).How these protocols work in detail is explained in [28,29].The proportion of the ACs undergoing SP-1 is P = 0.2, and the rest of them follows SP-2.By combining two SPs, peak reduction can be achieved.It can be observed that from 100 to 160 min, and from 400 to 600 min, the power demand of ACs based on the hybrid protocol combing two SPs has been tracking the reference output closely.It responded to the signal instantly and maintained the low power at the expected value for an hour (the natural TCL cycle time), but in the remaining time, the green curve gradually deviates from the reference value, and is slowly saturating close to the initial value.In contrast, at all times during the control period (100-350 min and 400-600 min), the power demand curve based on the proposed model and control strategy (red) is tightly maintained.
Although there are many advantages, the SP strategy is an open loop approach, and cannot compensate for the output error.On the contrary, the proposed closed loop control can adjust the control signal according to the results to obtain a higher precision.

Scenario 2: Simulation under Variable Ambient Temperature Conditions
In this scenario, we consider a more realistic situation, which means an ever-changing ambient temperature.We choose 5 August 2014, which was a hot summer day, and the total load was very high.The temperature sensors of the weather station deployed in the area collect practical ambient temperatures once every 10 min.We draw the ambient temperature profile of 24 h as shown in Figure 18.It is observed that from 12:00 to 15:00, the ambient temperature is high, and the maximum temperature is 36.3°C. Figure 19 represents 24 h power load of a downtown area in Nanjing, China, where there are 18,000 ACs of 10,000 residents operating on the hot day 5 August 2014.The real data for the load profile, which was collected every 15 min, was provided by State Grid Nanjing power supply company.The real-time temperature on the same day, which was collected every 10 min, was provided by the Jiangsu Provincial Meteorological Bureau.In Figure 19, the top blue curve is the 96 points load characteristic curve of the customers in the area on the maximum load day.We simulated the load values each second, by using the method of spline interpolation.Peak hours are from 10:30 a.m. to 16:00 p.m. and from 20:00 p.m. to 22:00 p.m.During these periods, the electricity demand outstrips supply, and a power shortfall will appear.As can also be seen, a low load period is from 4:00 a.m. to 8:00 a.m.The maximum load is 5.75 MW and the minimum load is 3.51 MW, and the load peak and off-peak difference is 38.96% of the maximum load.The bottom red curve in Figure 19 is for the baseline load of 18,000 ACs.The curve is obtained by simulating the two order linear time invariant transfer function model, which is based on the parameters in Table 1 and the ambient temperature curve in Figure 19.The total load characteristic curve minus the AC load curve based on the TF model is the green curve of the non-air conditioning load, which is in the middle of the figure.
The population consists of 18,000 ACs.Assuming that at the initial moment, 44.6% of the ACs are on, while the rest of them are off, the whole load control scenario is divided into two stages.The first stage is from 12:00 p.m. to 16:00 p.m., when the grid is in a peak load period.We decide to cut the peak load of ACs to 2 MW.The load aggregator will raise the temperature set point of the ACs to achieve peak load shifting.The second stage takes place between 4:00 a.m. and 8:00 a.m., when the peak load period of the power grid is over, we decide to increase the load of ACs to 1.2 MW to maintain the power consumption of the trough period.The load aggregator will drop the temperature set point of the ACs to realize users' comfort recovery.
In this scenario, the two order transfer function model of air conditioning is no longer linear time invariant, and instead, it is changing with time.The model parameters vary with the outdoor temperature in the case of equivalent thermal resistance, thermal capacitance and other parameters remaining unchanged.There are two ways to solve this problem.This paper focuses on the load control of ACs to achieve peak load shifting.Since the peak load duration is not long, the first method takes the average of the ambient temperature at intervals of half an hour within the required period of peak load shifting.In less than half an hour, it can be assumed that the ambient temperature is constant.During the period when ACs are participating in demand response, the change of the ambient temperature is not dramatic, especially for the existence of the heat island effect in the city center area.Therefore another method calculates the parameters of the transfer function whenever the outdoor temperature changes 0.5 °C.If the temperature variation is within 0.5 °C, the transfer function of the model is believed to remain unchanged.
In this paper, we adopt the second method.According to the different ambient temperature, the two order transfer function parameters are shown in Table 4.The numbers in Table 4 are subsequently calculated by using Equations ( 21)- (27).The relevant AC parameters used in the calculation can be found in Table 2.In the process of calculation, the ACs' other parameters are as follows: .The period when ACs could be controlled by the load aggregator is from 4:00 a.m. to 8:00 a.m. and from 12:00 p.m. to 16:00 p.m. Thus, according to the real-time ambient temperature on the warmest day, 5 August 2014, the research scope is identified as (28 °C, 32 °C) and (34 °C, 36 °C).
Figure 20 illustrates the effect of the AC load control based on the PSO tuning parameters of the PID controller for load shifting on a hot day when the ambient temperature continues to change.In Figure 20a, the upper blue curve represents the aggregated power 18,000 ACs without any control, the middle red dashed curve is the desired AC load reference value for load shifting, and the solid black line at the bottom represents the actual output power of AC load based on PSO tuning parameters of the PID controller.Figure 20b,c is obtained by zooming in on the parts circled by a dotted curve in Figure 20a.They show the effects of increasing the power output from 4:00 a.m. to 8:00 a.m. and reducing the power output from 12:00 p.m. to 16:00 p.m, respectively.Figure 20d shows the probability density distribution of the tracking error between the actual output and the reference value.Based on the discussions above, one can draw the conclusion that through the PID controller based on the PSO algorithm, the output power of 18,000 ACs can be controlled by adjusting the temperature set point offset.Throughout the control period, the output power of the air conditionings has been following the reference value closely, and the tracking error is small.During peak hours, the load aggregator can reduce the power of ACs by up to 6 MW, and maintain the total load of n ACs at about 20 MW.To some extent, the power shortage has been eased.In peak-off hours, the load aggregator can increase the power of ACs by 4 MW to make the load curve smoother, and maintain the total load of n ACs at about 12 MW.
Figure 21 shows the variation of ACs' temperature set point offset, which is also the control signal of the two order transfer function model for the air conditioning group.The range of air conditioning set point temperature is [−0.4 °C, 0.95 °C]. Figure 22 shows the variation of the indoor temperature of 30 random ACs.The hysteresis width of the room temperature, which varies with the air conditioning temperature set point u(t), has remained unchanged H = 1 °C.The figure illustrates such a process: the ACs' temperature set point can be controlled by the load aggregator who sets a higher value to reduce the load as well as setting a lower value to increase the load.The variation range of the indoor temperature is  At 16:00, when the load control period is over, the temperature set point offset does not return to 0. Instead, it continues changing (in Figure 21) and the room temperature does not directly return to the initial value [19.5 °C, 20.5 °C] (in Figure 22).The reason for this phenomenon is that, within the short time for the air-conditioning chillers to restart after shutdown, the power demand is usually significantly higher.This is the pattern of the air conditioning which makes the load rebound before and after the load control event.In order to bring the load curve back to the state before controlling, such as the curve after 16:00 p.m. in Figure 20a, a continuous control is required for some time after the load management.The novelty provided in this section is as follows: in the variable ambient temperature scenario, the two order transfer function model of air conditioning is no longer linear time invariant, and instead it changes over time.The model parameters vary with the outdoor temperature.This paper puts forward two suggestions to solve this problem, and we adopt the second one.

Conclusions
This paper presented the design of a closed-loop control strategy for air conditioning loads to participate in load control.A core contribution in this paper is to present the aggregated power demand of the AC population by describing how the proportion of operating ACs varies over time, in response to a step shift in the temperature set point of all ACs.Assuming that the thermal capacitances are distributed log-normally among ACs, a formula for the transient response to a temperature set point offset was derived.The temperature and power of ACs experienced an underdamped oscillation after the temperature set point had risen by 0.5 °C.By observing this phenomenon, we can use a simplified second-order dynamic model to approximate the air conditioning group.As far as we know, this is the first work to characterize a mathematical approximation for the period of such response.Based on the above work, a simplified control model based on the two order linear time invariant transfer function was derived.
Another contribution of this paper was the application of the closed-loop control and feedback strategy to provide a new way for flexible loads like air conditioning to be involved in demand response.In the feedback control strategy, the error between the expected demand reference value and the actual power of the AC group were used as an input variable for the controller.The output of the controller was the offset of the temperature set point, which was also the control signal of the AC group.
A simulation model for controlling the offset of the temperature set point was developed to verify that the proposed transfer function model and control strategy can closely track the reference peak load shifting curves.Two scenarios were selected for the simulation, including the constant ambient temperature and the variable one.The simulation results demonstrated the effectiveness of the controller design.The control model could lower electricity demand during peak hours effectively, fill the gap between peak and off-peak loads, and balance the supply and demand of electricity.The proposed closed-loop control strategy provided a new way for flexible loads to participate in demand response. .According to Section 3.1, the duty cycle is 0.5, so the ratio between the distance of two successive peaks that we call Proposing z1 = 0.9 and applying Equation (A4), we have: Besides, from the final value theorem in [30], we obtain: Assuming that the particle swarm contains i particles, the information of the particle i ( {1, 2,..., }) i l  can be represented by the D-dimensional vector which indicates the number of parameters that need to be optimized.i x = 1 2 ( , ,..., ) x represents the space position, and i v = 1 2 ( , ,..., ) represents the speed.After the two optimal solutions , i best p and , i best g are obtained, the particle swarm updates the velocity and position according to Equations (B1) and (B2): In the above formulas, ( ) i v t is the speed of particle i in the D-dimensional space at the time t, and ( ) x t is the position of the particle i in the D-dimensional space at the time t., ( ) i best p t represents the optimal solution of a single particle i itself, and , ( ) i best g t represents the optimal solution of the whole group. 1 c and 2 c are the acceleration factors and their general values are between (0, 2). 1 r , 2 r are random functions, scaling values range from 0 to 1.  is a non-negative weight, which affects the overall optimization ability.Figure 10 shows the flowchart of the PSO algorithm tuning PID controller parameters.The complete optimization procedure is as follows: Step 1: Initialize the PSO's population size, maximum number of iterations, the learning factor, as well as the initial position and velocity of particles.
Step 2: Choose the integrated time and absolute error (ITAE) as the fitness function.Calculate the fitness of each particle, find out the best individual in the initial particle swarm optimization, and initialize it as the best particle , ( ) i best g t in the group.Besides, the fitness of the particles themselves are taken as the initial value of the particle's individual optimal , ( ) i best p t .
Step 3: If the current fitness is better, take the current position of the particle as the best position , ( ) If the fitness of the particle's optimal position is higher than that of the swarm's optimal position, then take the optimal position of the particle as the best position , ( ) i best g t of the population.
Step 4: Adjust the speed and position of each particle according to the Equations (B1) and (B2).
Step 5: The terminating condition is that a predetermined maximum number of iterations or the lower limit of fitness has been reached.Otherwise, go to Step 3.
Step 6: Output the global optimal particles which are best parameters of the PID controller.

Figure 1 .
Figure 1.Periodical properties of the temperature inside the house and the power of air conditioners (ACs).
,b indicates the temperature distribution before the step response.Figure 2c,d illustrates the temperature change after the step response.Therefore, the temperature boundary   ,

Figure 2 .
Figure 2. The temperature distribution of ACs before and after a step change of 0.5 °C of the temperature set point.

0 Figure 3 .
Figure 3.The state distribution of ACs instantly after the temperature set point step change.

Figure 4 .
Figure 4.The average temperature and power for the population composed of 10,000 ACs before and after the step signal of temperature set point, respectively under conditions of rel  = 0.02, 0.05, 0.2 and 0.5.rel  is the standard deviation of log-normal distributions as a fraction of the mean value for R, C and P. (a) rel  =0.02;(b) rel  =0.05; (c) rel  =0.2;(d) rel  =0.5.

Figure 5 .
Figure 5.The aggregated demand response of a 0.5 °C step based on the second-order linear time invariant (LTI) model.

Figure 6 .
Figure 6.The dynamic performance of the transfer function.By analyzing the stability of the transfer function p G , we get the Bode diagram.As shown in Figure 7, the system amplitude stability margin p G and phase angle crossover frequency cg W are defaults, and the phase margin angle m P = −107.9061,gain crossover frequency cp  = 0.0227.Therefore, the closed-loop system with stability parameters marked in the figure is stable.

Figure 7 .
Figure 7. Bode diagram of the system transfer function p G .

Figure 8 .
Figure 8.(a) The variation on the average temperature inside the room and (b) the variation on the aggregated demand of the AC group when the ACs' temperature set point has a step rise of 0.5 °C.The two order linear time invariant transfer function model is compared with the ETP model which contains 50, 1000 and 10,000 ACs.

Figure 9 .
Figure 9.A load control strategy for the AC group based on the particle swarm optimization (PSO) algorithm for Proportional-Integral-Differential (PID) parameter tuning.

Figure 15 .
Figure 15.(a), (b), (c) A comparison between the desired reference value of the aggregated demand of a population of ACs for load shifting, the power output of the population based on conventional PID controller, and the power output of the population based on PSO tuning PID parameters, when the outdoor temperature is constantly 32 °C.(d) The aggregated demand of n ACs in the extreme case of a 100% load reduction for 30 min.

Figure 16 .
Figure 16.Variation of the indoor temperature of 30 random ACs in the population.

Figure 19 .
Figure 19.24 h power load of a downtown area of Nanjing, China where there are 18,000 ACs of 10,000 residents operating on the hot day 5 August 2014.

Figure 20 .
Figure 20.(a), (b), (c) A comparison between the aggregated power of 18,000 ACs without any control, the desired AC load reference value for load shifting, and the actual output power of AC load based on PSO tuning parameters of the PID controller.(d) Probability density distribution of the tracking error between the actual output power and the reference value.The model is simulated for a hot day, and the ambient temperature continues to change.

Figure 21 .Figure 22 .
Figure 21.The variation of ACs' temperature set point offset.

,
above are parameters of the denominator of the transfer function.Next we calculated parameters of the numerator.D(t) is the expected value of the expression the system is in steady state.In which on T and off T depend on the reference temperature Therefore, from the initial value theorem that is mentioned in[30], we get: Laplace transform of a derivative and the initial value theorem, we obtain:

Figure B1 .
Figure B1.The flowchart of PSO algorithm tuning PID controller parameters.

Table 1 .
Parameters of two order linear time invariant transfer function.

Table 2 .
Parameters of the air conditioners.
n  

Table 4 .
Parameters of two order transfer function model for ACs. )