The proposed multi-objective function comprises one economic, two technical, and one user-convenience functions, as introduced and discussed below:

To consider the effect of distributed energy resources such as wind and solar units, all OF _j are calculated in every possible scenario (S _i ) and then combined in the form of O F j = ∑ S 1 S n P S i ⋅ O F j S i 100 (7) where OF _j ^S _i and P(s _i ) are respectively the values of OF _j in scenario S _i and the probability of scenario S _i .

The constraints for the operation of the thermal power plants, load-shedding limits, as well as operation of ESSs and distributed energy resources should be considered for the successful implementation of the DR, as introduced and discussed below:

The power balance constraints of the power grid include P i j , t = Re S i j , t = V i , t 2 Z i j cos θ i j − V i , t V j , t Z i j cos δ i , t − δ j , t + θ i j (24) Q i j , t = Im S i j , t = V i , t 2 Z i j sin θ i j   − V i , t V j , t Z i j sin δ i , t − δ j , t + θ i j − b V i , t 2 2 (25) P i , t g + P i , t L S + P i , t W + P i , t S + − P i , t c h + P i , t d c h − P i , t L = ∑ j ∈ Ω l i P i j , t : λ i , t P (26) Q i , t g + Q i , t L S + Q i , t W − Q i , t L = ∑ j ∈ Ω l i Q i j , t : λ i , t q (27) in which (24) and (25) describe the relationship for respectively active and reactive power in the network while (26) and (27) outline the active and reactive power balance constraints at time t. Also, P _i,t ^L is the power demand and Q _i,t ^W is the reactive power generated by the wind power plants on bus i at time t.

The DR constraints are formulated as follows: P i , t L = P i , t L 0 + D R i , t ;   ∀ i (28) ∑ t = 1 24 D R i , t = 1 ;   ∀ i (29) − α ⋅ P i , t L 0 ≤ D R i , t ≤ α ⋅ P i , t L 0 ;   ∀ i (30) tan θ i , t = Q i , t L 0 P i , t L 0 (31) Q i , t L = ⁡ tan θ i , t P i , t L (32)

Constraint (28) specifies that the electrical power consumption at bus i and time t in the presence of the DR scheme (P _i,t ^L ) which is obtained from the sum of the base load (P _i,t ^L0 ) and decision variable of the DR (DR _i,t ). When DR _i,t is zero, no DR program occurs, however for DR _i,t >0, the active and reactive power demand of bus i will increase at hour t, resulting in valley filling. Similarly, the active and reactive power demands of bus i decrease at hour t when DR _i,t <0; hence, contributing to peak shaving. In accordance with (29), the net change in the demand across all buses is zero. The upper and lower bounds for the DR values are determined by (30). By maintaining a constant power factor before and after DR, as described by (31), the coordinated adjustment of the reactive power in the DR scenario is expressed by (32).

The relationships between the HVAC units can be simplified by omitting the coefficient related to the heating system and expressing the coefficient and value of the cooling system in an aggregated form, as described in (33). This simplification aims to reduce the volume of calculations and assumes that the calculations are performed on a hot day of the year. In the non-flexible mode, the ambient temperature, room temperature, and air conditioning system capacity are considered; however, in the flexible mode, apart from these parameters, the value of k ₃ . P _AC (i,t) is determined based on the load chart specifications and the specified ranges of the peak and valley loads. Constraint (34) specifies that intelligent control is employed to ensure that at low loads, the electricity consumption is higher and the temperature is maintained close to the minimum allowed temperature range. The goal is to minimize electricity consumption and adjust the temperature to the maximum allowable temperature range. T i n i , t = k 1 . T i n i , t − 1 + k 2 . T a t − 1   − k 3 . P A C i , t . U i , t (33) U i , t = 1 T i n i , t − 1   > T Max   1   T a t > T a t − 1   1 k 1 . T i n i , t − 1 + k 2 . T a t − 1   > T Max   0 e l s e   (34)

In this study, the method borrowed from (Soroudi, 2017), which was used in the DC-OPF program. The allocation of ESSs and their sizes minimized the operational cost. After demand prediction and HVAC aggregator calculation in the AC-OPF platform, the formulation size and number of ESS are calculated. Considering the power quality issues, the uncertainties associated with renewables are addressed when optimizing ESS placement. The state of charge (SOC _i,t) represents the charging capacity of each ESS, depending on the number of ESS units installed on-site (N _i ^ESS ) and their minimum and maximum capacities (SOC _i,min and SOC _i,max ), as described by (35). Constraint (36) specifies that the total number of installed ESS units should not exceed the specified maximum limit (N _max ^ESS). S O C i , ⁡ min × N i E S S ≤ S O C i , t ≤ S O C i , ⁡ max × N i E S S (35) ∑ N i E S S ≤ N max E S S (36)

The uncertainties mainly arise from random variations in the input data, prediction errors, network failures, and intermittent behavior of RESs. In such an uncertain environment, deterministic power flow (DPF) is infeasible and cannot precisely reveal the state of the system. Therefore, for uncertainty assessment, several uncertainty modeling techniques exist, such as probabilistic approaches, possibilistic approaches, hybrid probabilistic-possibilistic approaches, information gap decision theory (IGDT), interval analysis, and robust optimization (Singh et al., 2022). One of the possibilistic methods is scenario-based. Although the scenario-based method has disadvantages such as sensitivity to selected scenarios and complexity in the number of scenarios, advantages such as flexibility and the possibility of performing sensitivity analysis have led to the use of this method. The three scenarios of normal, high, and low production for each solar and wind unit are shown in Figure 6. Each of these three scenarios was assigned a probability of occurrence in such a way that for the three scenarios of normal, high, and low production in the wind system, the probability of occurrence is 50, 25, and 25, and in the solar system, it is 60, 20%, and 20%, respectively. Owing to the existence of three scenarios for each of the sources, and without losing the generality of the problem, nine scenarios were obtained, and their probabilities are depicted in Figure 7. In each of the nine scenarios, the problem was solved and multiplied by the probability of its occurrence, and the system output was obtained from the sum of the results of the nine scenarios, taking into account the uncertainty of scattered production resources.

Accurate electricity demand prediction is vital for efficient network operation because errors can result in increased production costs and grid losses. Also seasonal ARIMA has better forecasting result for capturing the seasonality and short-term fluctuations in demand but in general, LSTM performed better for products with stable demand (Falatouri et al., 2022). In this study, the Modified LSTM method for load predictions, based on the deep neural network architecture referenced in (Zarei and Ghaffarzadeh, 2023), was utilized. The prediction data cover the period from August 19 to 19 September 2022, and were obtained from the New Zealand Electricity Company (Live load data, 2022) using the MATLAB software. Figure 8 demonstrates the load prediction using this method, with the achieved accuracy indicated by the root-mean-square error and mean absolute error values of 61.8 and 274.9, respectively.

Two flexible and non-flexible modes were considered for the aggregator to investigate the effect of the HVAC aggregator on DR. As mentioned earlier, this aggregator is included in Bus 8 with a capacity of 0.4 megawatts and a power factor of 0.8. In flexible mode, the aggregator adjusts its performance according to the predicted load chart. In addition to the ambient temperature inside and outside the building, it also checks the load chart and then makes a decision regarding the settings of the on and off times and the corresponding power. For a better comparison, in addition to the simultaneous effect of flexible aggregators and non-flexible, different percentages of DR were applied to other buses. In Figures 9–12, respectively, the objective functions of the total cost, total voltage variation, peak reconnection, and use convenience are displayed in two different modes of aggregator operation and in different DR steps. As can be seen in all the different DRs, the total system cost and total voltage variation in the flexible operating mode of the HVAC aggregator are reduced, and the amount of peak reduction also increases in this mode. The aggregator installed at a point does not have a noticeable effect on the fourth objective function, that is, the UC. To investigate the effect of the HVAC aggregator individually, in the state without DR, the demand of the entire network in two flexible and nonflexible modes is shown in Figure 13. Because its capacity is approximately 10% of the demand of the entire system, a noticeable reduction in the two peak loads and the transfer of its consumption time to periods of low and intermediate loads are observed.

According to the possible scenarios for RESs, the total cost of the system in different DRs from zero to 10% is calculated and shown in Figure 14. Therefore, in scenarios with the possibility of producing more energy through solar and wind energy, the cost of the entire system is reduced. In addition, the cost of the entire system decreased significantly with an increase in the DR. Scenario 9, characterized by maximum solar and wind production, exhibits the lowest cost at each DR step, whereas Scenario 5, with minimal renewable resource production, incurs the highest system cost.

In Figure 15, the voltages of all buses are shown under peak load and low load conditions, with and without consideration of the DR. As can be seen, the effect of applying the DR is to improve the voltage of all buses and approach the optimal voltage (1 pu), especially in low load and end buses of the distribution network.

The effect of levelization on the demand graph of all network loads by applying DR is presented in Figure 16. In the case without DR, the load curve had a desert valley and almost two afternoon and evening peaks. By applying only 10% DR, valley filling has been done to a considerable extent, and almost all the two peaks mentioned have also had a significant reduction. This important factor causes fewer negative effects on the smart distribution grid during peak times.

In this subsection, the load scheduling involving demand response at various generating stations, including conventional, wind, solar, DG, and ESS, is discussed. Additionally, the electrical load power, including electrical load and power loss, is presented both without DR and with 10% DR states. As depicted in Figure 17, the implementation of a 10% DR results in a reduction of conventional generation and a subsequent decrease in generation costs. Furthermore, with the 10% DR, power generation costs are lower at lower loads, leading to an increase in power generation by DGs compared to scenarios without demand response.

In order to address multi-objective optimization challenges, various methods like the Ɛ-Constraint and TOPSIS are utilized. Each method has its own set of advantages and drawbacks. Figure 18 illustrates the values of four objective functions: total cost (TC), total voltage variation (TVV) for power quality enhancement, peak reduction (PR), and user convenience (UC), ranging from 0 to 30 percent of DR. Utilizing the TOPSIS method identified the optimal point in DR at 30%, albeit resulting in a significant decrease in user convenience from 100% to nearly zero. Given the opposing trends of the other three functions in comparison with UC, even with low weighting factors, the final solution remains unaffected. To address this issue, a thorough analytical approach to the objective function responses is imperative. The cost function exhibited a fully decreasing trend, with a steeper decline at lower DR values and a more gradual decrease at higher DR values. The TVV objective function displays a slight decrease with increasing DR, with a reduction range of less than 0.7% between 23.7 and 23. Because a substantial reduction in TVV comes at a high cost, a balance must be struck. Peak reduction, as the third objective, sees a rapid increase in peak value at the start of DR, reaching a local peak of over 10% before 15% DR, followed by minimal changes thereafter, and maintaining peak reduction below 12% up to DR 30%. Despite increasing the DR in certain regions, a decrease in the PR may still occur. Finally, the UC objective function focuses on user convenience, which diminishes gradually with DR, drops sharply after 15% DR, and approaches zero at 30% DR. Considering the explanations provided, the importance of subscriber consumption patterns, and the societal impact of economic decisions, an analytical approach suggests that DR values below 15% are optimal for altering subscriber consumption patterns. During this phase, when UC falls below 20%, significant cost reductions across the system and optimal peak efficiency are achieved at a DR of 14%.

The simultaneous solution of AC-OPF and DR problems in a smart distribution network was achieved by applying the effects of a small-scale load aggregator owing to its control and feedback capabilities. Although, in theory, there is no limit to the participation of more electricity subscribers in demand response, in practice, a large change in the time and amount of their consumption will cause dissatisfaction. In practice, choosing the appropriate percentage of DR by considering the behavioral habits of energy consumption of electricity subscribers in the framework of user convenience makes the network operator’s planning on this value closer to reality and reduces the difficulty of adjusting the demand-production balance. The participation of the local aggregator of the HVAC loads to comply with the optimal temperature range of the residents was also emphasized.