DC Nanogrids for Integration of Demand Response and Electric Vehicle Charging Infrastructures: Appraisal, Optimal Scheduling and Analysis

: With the development of electronic infrastructures and communication technologies and protocols, electric grids have evolved towards the concept of Smart Grids, which enable the communication of the different agents involved in their operation, thus notably increasing their efﬁciency. In this context, microgrids and nanogrids have emerged as invaluable frameworks for optimal integration of renewable sources, electric mobility, energy storage facilities and demand response programs. This paper discusses a DC isolated nanogrid layout for the integration of renewable generators, battery energy storage, demand response activities and electric vehicle charging infrastructures. Moreover, a stochastic optimal scheduling tool is developed for the studied nanogrid, suitable for operators integrated into local service entities along with the energy retailer. A stochastic model is developed for fast charging stations in particular. A case study serves to validate the developed tool and analyze the economical and operational implications of demand response programs and charging infrastructures. Results evidence the importance of demand response initiatives in the economic proﬁt of the retailer.


Introduction
The Smart Grid paradigm requires active participation and coordination of all the agents involved in the electric system [1,2]. With the advent of this concept, nanogrids (NGs) have emerged as an invaluable framework for the integration of renewable generation, electric vehicles (EVs), storage facilities and demand response (DR) programs [3][4][5][6]; as well as for electrical supplying of remote isolated areas [7]. In this context, the optimal coordination of the agents involved requires deploying advanced communication infrastructures and electronic interfaces on either AC, DC or hybrid networks [8][9][10]. Communication channels enable effective information exchange among the different assets that participate in the operation of the system. This novel paradigm allows to effectively coordinate the participation of the different agents, with the aim of pursuing more efficient management of the generation and storage facilities [11]; while consumers can also actively participate through DR initiatives, modifying their consumption profiles on the basis of price or incentive signals [12].
Because of the heterogeneity and conflicting interests of the different agents involved in the NG operation, their coordination is frequently addressed in a centralized fashion by the grid operator. This agent is usually integrated into the local service entity along entity structure along with the local retailer. This way, in contrast to most of the existing scheduling tools for similar grids which are devoted to minimizing the operational cost, the developed formulation aims at maximizing the benefits obtained by the retailer optimally coordinating on-site generators, storage facilities, flexible consumers and public charging infrastructures. In addition, the proposed model incorporates various uncertainties, which are modeled as stochastic parameters, including the charging station demand for which, in contrast to other approaches which model the stochastic EV charging demand using conventional probability functions, a specific stochastic model is developed. For the sake of clarity, the major contributions of this paper are summarized below: • Describing a DC layout for isolated NGs with DR programs and public EV charging stations.

•
Developing an optimal day-ahead scheduling model for the grid system described, which is suitable for a grid operator which is integrated into the local service entity structure along with the local retailer. • Developing a stochastic model for the different uncertainties involved in the NG operation, including renewable generation, local demand and EV charging profiles.

•
Analyzing the influence of DR programs and EV penetration in the NG operation, also focusing on the economic impact of such aspects.
In the rest of the paper, Section 2 describes the analyzed NG layout, the agents involved and their role in NG operation. Section 3 presents the developed optimal scheduling model for the NG described in the previous section. Section 4 develops a stochastic framework for considering uncertainties in renewable generation, demand and EV charging profiles. Section 5 presents various case studies. The main conclusions are duly drawn in Section 6.

Description of the NG under Study
As commented, this paper is focused on isolated DC NGs. This kind of system involves the participation of various agents. The more notable interactions among the agents involved are depicted in Figure 1, while their roles are explained below: • NG operator: this agent is integrated into an upscale structure called a local service entity. It is responsible for operating the grid in an optimal way, ensuring the supplying quality and reliability. To this end, this agent daily performs a day-ahead optimal scheduling plan, by which the different on-site assets are coordinated with DR premises, such as those enabled by flexible consumers and public charging infrastructures. As a result, commitment signals, power set-points and DR information are sent to generators, storage systems and flexible consumers to address the resulted scheduling plan. • Retailer: this agent provides fuel for conventional generators and is responsible for the monetary expenditures derived from generator operation (operation and maintenance). On the other hand, it receives monetary incomes from consumers through energy tariffs and public charging infrastructures, of which the local service entity is the owner. • Generators and storage facilities: they are on-site assets owned by the local service entity. They may be formed by conventional generators, such as diesel engine generators (DEGs); renewable generators, such as PV panels or wind turbines (WTs); and storage facilities like BES. • Consumers: residential demand and public charging infrastructures are considered as consumers in this paper. The residential demand comprises inelastic and flexible consumption. While the first one does not respond to price or incentive signals from the NG operator, the second one may be scheduled in order to increase the efficiency of the system or ensure its reliability. On the other hand, public charging stations are owned by the local service entity. They provide adequate charging infrastructures to privately owned EVs, obtaining a monetary counterpart which is received by the retailer. To effectively address the optimal coordination of the different participants, both energy and communication channels have to be enabled, which are illustrated in Figure 2. Each controllable agent is connected through electronic interfaces, which allow an effective energy flow control from/to the different elements. This way, PV arrays, charging points and BES are directly connected to the DC bus through DC-DC converters. On the other hand, DEG and WTs require a rectified stage, which converts the generated AC waveform to DC. It is assumed that residential consumers mostly comprise DC loads; therefore, they can be connected directly to the DC bus without the necessity of an interface. In contrast, flexible consumers require a control mechanism that receives the commitment signals emitted from the NG operator. Each day, the operator receives forecasted data for weather and demand, which allow the performance of an optimal scheduling plan for the grid. This plan is transferred to the conventional generators and storage systems in the form of commitment and set-point signals, while renewable sources receive power setpoints. DR programs are enabled through direct communication with flexible consumers. Thereby, the public charging station may be controlled through power set-points, so that EV demand may be totally satisfied or not, while flexible residential consumers can be actually scheduled, shifting their consumption to the more convenient hours of the day.

Optimal Scheduling Model for the NG under Study
In this section, the mathematical formulation for the optimal scheduling of the NG described in Section 2 is developed. The optimization problem is formulated as Mixed-Integer-Lineal programming, which can be solved using conventional solvers.

Assumptions
It has been widely reported that electronic devices can delay the response of the device to which it is interconnected by several seconds [37]. This aspect was not considered in the developed formulation because it was assumed that the time step scheduling is considerably larger than the delay introduced by electronic devices. This way, this aspect did not impact the scheduling result yielded by the developed tool.

Objective Function
The optimal scheduling problem is formulated from the retailer point of view, which aims at maximizing its own profit. To this end, this agent obtains monetary incomes from serving energy in flexible and non-flexible consumers as well as charging processes in EV stations. On the other hand, this agent incurs monetary expenditures due to fuel consumption of DEG and operation-derived costs of the other generation and storage assets. Keeping this in mind, the objective formulation can be formulated as follows: where R stands for the representative scenario set; ω r is the probability of the r th scenario; ∆τ is the time step; T is the set of time intervals; λ is the energy cost (energy tariffs); d is the expected demand; the superscript EV makes mention to electric vehicles; Q is the set of shiftable consumers; u is a binary variable that indicates the commitment status; p stands for power; ρ is the operation and maintenance costs; and the superscripts PV, WT, BES and DEG make mention to solar panels, wind generators, battery storage and diesel engine, respectively. The objective formulation is formulated considering a stochastic model for uncertainties so that the expression (1) considers the probability of occurrence of each scenario. The first term of the incomes stands for the energy consumption of inelastic residential consumers, the second term counts the profit from EV charging and the last element represents the payments of flexible consumers.
The monetary expenditures mainly comprise operation and maintenance costs of generation assets. For renewable-based generators, these costs are proportional to the energy generated [7], whereas, in the case of BES, they are quadratic functions of the energy exchanged with the grid [38]. The last element of the expenditures stands for the fuel cost of DEG, which can be calculated as a quadratic function of the energy generated, as follows [39]: where a DEG , b DEG and c DEG are predefined coefficients. In order to keep the lineal feature of the developed problem, quadratic functions can be efficiently linearized using piecewise representations [15] (see Appendix A). (3) and (4), respectively.

DEG units can be modeled as other dispatchable generators by their upper and lower-rated values and ramping constraints, as said the Equations
where p DEG and p DEG are the minimum and maximum dispatchable powers of the DEG and R is the ramping rate.

PV Modelling
PV generation is influenced by uncertain parameters such as solar irradiation and ambient temperature. Assuming these parameters are known or sufficiently well-predicted, the available power generation can be calculated using an available PV array model. In this paper, the model considered in [7,14] is used, which determines the PV generation capacity as a function of the solar irradiation and ambient temperature, as follows: where p PV is the rating power of the PV array; ϑ is the solar irradiance; θ is the ambient temperature and η PV is the efficiency of PV panels.
As pointed out in [7,14], model (5) does not take into account the installed power of the PV generators. In other words, the second term in (5) may be higher than 1, which yields a value higher than the rated power of PV panels. In practice, inverters impose limitations in the PV power generation to avoid failure or fast degradation of components. Traditionally, the equipment allows exceeding the nominal power by 10% [40], assuming that these overvalues may be eventually assumed without excessively deteriorating the components of the PV array. With these premises, the PV array model is completed with (6).
It is worth noting that the developed model considers the PV generation as a variable of the problem. It means that the NG operator can schedule its production on the basis of forecasted data. In this sense, the scheduling result may occasionally set the PV generation below the available power calculated by (5). In such a case, the surplus energy is assumed to be dissipated in dummy loads [41,42].

WT Modelling
As in the case of the PV panels, the power given by a WT depends on the available wind speed. Once this parameter is known or approximated, the available wind generation can be estimated on the basis of the wind speed-power curve of WTs [7]. A typical graph for this kind of curve is shown in Figure 3. It can be distinguished into three clear zones. For a wind speed below γ WT , the turbine is not able to extract energy from the wind and therefore, the power produced is zero. After surpassing this threshold, the power given by the turbine increases until reaching the rated power of the device at γ WT, * , which is a characteristic value of each turbine model. Beyond this point, the turbine is able to give its rated power until a maximum threshold γ WT , beyond which the turbine is braked to avoid breakdowns. The wind speed-power curve of a WT is typically represented by a set of equations as shown in (7), while the expression (8) takes into account the efficiency of the coupling rectifier.
where γ is the wind speed; γ WT and γ WT are the minimum and maximum wind speed values for which the WT can work; γ WT, * is the nominal wind speed of WT; α WT and β WT are predefined coefficients; p WT is the rated power of WT; and η WT is the efficiency of WT.

BES Modelling
The power that a BES can exchange with the system is typically upper bounded, as said the constraint (9), whose nominal values are determined by the nominal capacity and the energy to power ratio [43]. It is realistic to assume that charging and discharging processes are actually complementary, which is ensured by imposing the constraint (10). Equation (11) models the state of charge of the batteries, which is limited by their nominal capacity and depth of discharge settings [40], as said the Equation (12).
where ε BES is the energy stored in batteries; η BES is the efficiency of batteries; ε BES is the nominal capacity of the storage system; DOD BES is the depth of discharge; and the superscripts BES, ch and BES, dch stand for the charging and discharging processes of batteries, respectively. Equation (11) models the energy stored in batteries as a function of the state of charge at t − 1 and the total energy exchanged with the grid. However, as pointed out in other studies [7,40], this model has to be completed by setting the state of charge at the beginning of the time horizon. As customary (e.g., see [41]), it is assumed in this work that the BES is fully charged at the beginning of the time horizon. In addition, to keep the model coherent, the state of charge of the BES is forced to be equal to the nominal capacity at the end of the time horizon. These working conditions are ensured by imposing the constraint (13).

Shiftable Consumers Modeling
As commented previously, it is considered that a part of the local demand may respond to commitment signals emitted by the NG operator. This kind of responsible load can therefore adjust their scheduling programs to achieve a more efficient operation of the grid. As in other related problems [14], it is assumed that these consumers have to be operated continuously, which means, once their operating cycle has started, it cannot be interrupted until completing their duty cycle, as modeled the Equation (14).
where on and o f f are binary variables that indicate the activation/deactivation of a shiftable consumer. It is assumed that these consumers have some preferences about when they prefer to be scheduled, in this sense, they have to be operated within allowable time windows, the scheduling tool being free to schedule their consumption in those time ranges. The constraint (15) ensures that shiftable consumers complete their duty cycles within the considered time windows. The model has to be completed by Equations (16) and (17), which ensure that responsible loads cannot be scheduled out of their time windows and are activated once over the time horizon, respectively. ∑ t∈Ψ q u q r|t = ϕ q ; ∀r ∈ R ∧ q ∈ Q (15) ∑ t/ ∈Ψ q u q r|t = 0; ∀r ∈ R ∧ q ∈ Q (16) ∑ t∈Ψ q on q r|t = 1; ∀r ∈ R ∧ q ∈ Q (17) where ϕ is the total number of time slots that a shiftable consumer must be connected throughout a day. On the other hand, and with the aim of completing the analysis, it has been also considered the case in which the shiftable loads could not be operated in a flexible manner. In such a case, it is assumed that these consumers desire to start their operation at the first slot of the allowable time window, which is ensured by the model (18)- (20).

Public EV Charging Station Modeling
The developed model assumes the public charging station can be operated in a flexible manner, however, in contrast to the flexible consumers considered in the previous section, the charging points response in this case to set-point signals rather than commitment orders, as said in Equation (21).
By the model (21), the expected EV demand may be only partially satisfied. This way, the operator has the capability to decide if it is necessary to not fully cover or even ignore some charging events, whether this supposes an actual benefit for the grid operation.

NG Balance
The energetic balance is ensured each time step by the constraint (22), which establishes the generation-consumption equilibrium. In this case, the non-served load was not contemplated, thus having to satisfy the inelastic expected demand entirely. Therefore, DR is in this case enabled by shiftable consumers and the public charging infrastructure.

Optimization Problem
The optimal scheduling model developed above is formulated from the retailer point of view, which aims at maximizing their profits while the operator ensures the quality of supplying. In this sense, the model is formulated as a maximization problem in which (1) supposes the objective function. Therefore, the result of the problem searches for the optimal scheduling of the different controllable assets that maximizes the monetary incomes of the retailer. For the sake of analyzing, a case was also considered in which the shiftable consumers cannot be operated in a flexible way, minimizing the DR capability of the system. Both problems are stated below.

Uncertainties Modeling
The optimal scheduling model developed in Section 3 contemplates various uncertain parameters, which are reported in Table 1. For managing such parameters, a stochastic framework has been proposed. Stochastic modeling is based on representing the uncertain parameters by means of probability distributions, from which one can generate a number of scenarios [44]. The scenarios generated constitute the scenario-space (S), which has to be formed by many members since a large number of scenarios have to be considered to properly catch the stochastic essence of the parameter (normally 1000 [45]). Because of the large size of the scenario space, the optimization model may be intractable in practice. Indeed, one should note that all the variables and constraints would have size S × T . To overcome such issues, some papers have proposed to use scenario reduction approaches [14,45]. By these techniques, the original space is completely described by just a set of few representative members, which are considered a sufficiently descriptive image of the original space. Among the available space reduction approaches, clustering techniques are quite popular because of their efficiency, adaptability accuracy and simplicity [46]. Specifically, the k-medoids method has been used in this paper because of its overall good performance [46]. The k-medoids method is inspired by the optimal plant location problem. In such a problem, a set of plants must be placed in a set of cities, so that the total distance from the plants in the cities is minimized. This approach thus grouping the set of scenarios into clusters, so that each cluster is fully described by just a member within it called medoid. This way, the scenario space is reduced to the so-called representative scenario-space (R), which is a notably smaller scale. The k-medoid method presents a degree of freedom, namely the total number of clusters to be created. To set this parameter, helpful indicators such as the total sum of distances and the Davies-Bouldin index along the methodology described in [47] have been used. The flowchart in Figure 4 summarizes the steps necessary for solving the optimal scheduling problem developed in Section 3 with stochastic modeling of uncertainties. One of the most interesting features of the k-medoids method is the possibility of easily calculating the probability of occurrence of each scenario, as follows: where S is the scenario space; and Ω r represents the cluster of the r th representative scenario.  It remains to be explained how the different uncertain parameters involved are modeled on the basis of probability functions, which is explained in the following sections.

Predictable Parameters
Some uncertain parameters can be day-ahead forecasted with sufficient accuracy [39]. Examples of these profiles are weather parameters or inflexible demand. For these kind of uncertainties, the forecasted data can be subjected to errors derived from the forecast technique or unpredictable events. Such errors can be modeled as normal distributions [48], taking the forecasting parameter as mean and a pre-set standard deviation. Using this distribution function, the scenario space can be generated for these parameters. For the model developed here in particular, the inflexible demand, solar irradiation, temperature and wind speed are assumed to be managed by this approach [38,49,50].

EV Demand Modeling
In contrast to the other uncertain parameters, the EV demand, especially at lowaggregated levels, is hardly predictable due to the generally random behaviors of drivers. In this sense, some references have considered specific probability distributions for EV demand [16,51,52]. For the particular case of the NG under study, it is important to note that at this scale the majority of consumers are domestic. In this sense, it was considered that the public charging station is devoted to light vehicles (cars and motorbikes).
In this paper, a stochastic model for the considered public charging stations has been developed. The developed model is based on the vehicle trip distribution (F), which is shown in Figure 5 for a typical weekday in US [16]. This distribution indicates the probability of an EV is driving on road at each hour of the day. It is realistic to assume that the probability of charging events will increase with the probability of a vehicle driving on the road [16]. Thus, the probability of charging events each hour of the day can be adjusted to the probability distribution shown in Figure 5. On the other hand, it is necessary to model the number of charging events that may occur in a day. It is assumed that this parameter can be either predicted based on historical data or directly estimated from for example market studies. In this sense, the number of daily charging events (N V ) is modeled with a normal distribution, where the mean µ stands for the daily expected number of charging events and σ models the standard deviation, as follows: N V s = round(rand(N(µ, σ))); ∀s ∈ S, ≥ 0 (26) where round(·) is a function that yields the nearest integer; rand(·) is a function that yields a random number based on a probability function; and N(µ, σ) is a normal distribution with mean µ and standard deviation σ.
Once the total number of charging events was determined, it can be used as the entry of the vehicle trip distribution, in order to determine the arrival time of each vehicle, as follows: T s = round(rand(F)); ∀s ∈ S, T i|s = T 1|s , T 2|s , . . . , T N Vs |s (27) Now, the vector T can be used to construct the vector c, whose i th element is equal to 1 if a vehicle arrives at the charging station at the i th time instant and 0 otherwise, using the following simple loop: initialize c s = c s|1 , c s|2 , . . . , c s|T = [0, 0, . . . , 0] ∈ B T for i = 1 : size(T s ) do c s (T s (i)) == 1 end do ; ∀s ∈ S (28) Now, the vector c can be used for modeling the EV demand, as follows: d EV s|t = c t P EV Rated ; ∀s ∈ S ∧ t ∈ T (29) where P EV Rated is the rated power of the charging process. The model (29) considers two plausible assumptions:

•
The EV demand is considered constant during all the charging processes which, as indicated in [15], is a quite realistic assumption since only marginal variations with respect to the rating values are observed during the short time of the charging event.

•
The EV charging process is completed within a unique time slot, which is quite realistic, assuming fast charging processes, which can be completed in only 15-30 min [15].
Nevertheless, it is worth mentioning that the developed stochastic model for EV demand presents a modular structure that allows it to be adapted to any charging mode or more comprehensive charging profiles. Nevertheless, it has been assumed that the developed model is quite realistic and useful in practice, lying out of the scope of this work further analyzing other more elaborated models. For the sake of example, Figure 6 plots some EV demand profiles constructed with the developed model considering various numbers of expected charging events. In this figure, fast charging mode was considered, with 55 kW rated power [15].

Case Study
This section presents a case study that is devoted to validating the developed stochastic optimal scheduling formulation, for the DC NG layout described in Section 2. The optimization model was coded in Matlab R2019a and solved using Gurobi [53] with 30-min resolution.

Data
For creating the scenarios for uncertain parameters, several public databases were used. Figure 7 shows the forecasted profiles for the solar irradiance, temperature, local demand and wind speed along the scenarios generated using the methodology described in Section 4. The considered weather profiles were taken from Ref.
[54] and correspond with real data observed at the Virgin Islands on 3 May 2016. The forecasted profile corresponding to local demand was adapted from real demand measurements at La Palma Island at the same date, which is available in [55]. The data corresponding to the other parameters involved are reported in Table 2, while the data regarding shiftable consumers are collected in Table 3. For the EV charging the fast mode was considered, with a rated power of 55 kW and a price of 1.5 $/kWh [15]. Finally, Figure 8 shows the Time-of-Use tariff used for non-flexible consumers, which is based on typical dynamic tariffs [14].

Results
In this section, various results are presented and commented on. For further analysis, the operation of the NG under study is considered with and without flexibility in shiftable loads (models 1 and 2, respectively). Moreover, the model is tested for various numbers of expected charging events (µ in Equation (26)). Figure 9 shows the value of the objective function (retailer profit) for various cases. As observed, the flexibility provided by consumers has a direct impact on economical profitability. In fact, the monetary balance for the retailer resulted in negative values (monetary losses) if the shiftable consumers are operated in a rigid manner, while the benefits may achieve up to~400$ in the case of the flexible program. Logically, the profitability of the system grows with the number of expected charging events, since more monetary incomes are obtained from the public charging station, as seen in Figure 10.  As seen in Figure 9, the flexibility of shiftable consumers has a direct impact on the profitability of the system. This is due to the flexibility provided by consumers allows to further satisfy the expected EV demand. Indeed, as observed in Figure 11, the EV demand satisfaction decreases in the case of considering the model, in which demand flexibility is not enabled. In order to further analyze the behavior of the charging station, Figure 12 plots the expected and actual satisfied EV demand with model 2, as seen in this figure. Some expected charging events were not covered, while others were just partially satisfied.  In order to further analyze the role of the shiftable consumers, their scheduling results are plotted in Figure 13 considering the optimization models 1 and 2. As seen, in the case of flexible operation, both consumers are scheduled later, exploiting higher PV penetration and avoiding the use of DEG. As commented, the main advantage of operating the shiftable consumers in a flexible way is reducing the dependency on costly generators like DEG. As seen in Figure 14, the expected DEG generation with flexibility barely surpassed 400 kWh, while it was higher than 600 kWh in case of no flexibility. These results are traduced in lower fuel consumption and therefore lower operational cost. As seen in Figure 15, the expected fuel cost with model 1 did not surpass 400$, while the fuel cost was higher than 800$ in the case of the rigid operation of shiftable consumers. These results have also a direct impact on environmental concerns since, as commented in [40], total CO 2 emissions are proportional to diesel consumption.

Conclusions and Future Works
This paper has comprehensively analyzed a possible DC isolated nanogrid layout for integration of renewable sources, energy storage facilities, demand response initiatives and electric vehicle charging infrastructures. In this sense, the necessary electronic interfaces and communication channels have been described. Moreover, the relationship among the different agents involved in the grid operation has been discussed. An optimal scheduling tool for the analyzed nanogrid has been presented, along with a stochastic framework for uncertainties modeling. In this regard, a specific stochastic model for fast charging stations has been developed.
A case study has been presented and various results have been commented, with the aim of validating the developed formulation and analyzing the impact of demand response programs in the grid operation. The results obtained have served to evidence the importance of flexible consumption for boosting up the monetary incomes of the retailer and reducing diesel consumption and therefore CO 2 emissions. It has been also remarkable the role of public charging infrastructures in the profitability of the grid and how demand response programs affect their viability.
Future works should be focused on analyzing the considered grid layout in gridconnected mode, for which the uncertainty of energy pricing along with the possibility of selling energy to the upscale grid has to be considered.

Acknowledgments:
The icons used in the figures of this paper were developed by Freepik, ultimatearm, monkik and smalllikeart from www.flaticon.com (accessed on 16 September 2021).

Appendix A. Linearization of Quadratic Terms
In the developed optimization model, some quadratic terms have appeared in the objective function. To linearize these terms and, therefore preserve the lineal character of the problem, a piecewise representation of such variables has been considered. By this approach, the quadratic function ψ is characterized by its piecewise representation ψ, for which n ( x i ) points are taken, as follows: ψ = x i , ψ( x i ); ∀i ∈ {1, 2, . . . , n} (A1) Usually, the higher number of points considered, the more accurate the piecewise representation is however, the optimization problem becomes computationally more expensive. The piecewise representation of quadratic terms allows to replace them by the dummy variable z wherever they appear in the problem. This variable is calculated, as follows: where δ is a binary variable that ensures that only one segment of the piecewise curve is activated at once by imposing the following constraint: To ensure that only one segment of the piecewise curve is activated at once, the following constraint must be imposed: In (A2), the parameters L and K can be calculated as follows: Finally, it is important to note that (A2) involves the product of two variables (δ and x), which introduces further nonlinearities in the model. To avoid this issue, the product of variables can be replaced by the lineal variable w for which the following additional constraints are needed: where M is a large positive constant.