Multisources of Energy Contracting Strategy with an Ecofriendly Factor and Demand Uncertainties

This study presents a mathematical formulation to optimize contracting capacity strategies of multisources of energy for institutional or industrial consumers considering demand uncertainties. The objective consists of minimizing the total costs composed of the different types of energy contract capacity costs, penalty price, and an ecofriendly factor. The penalty price is charged on the demand of energy exceeding the total contract capacities. The ecofriendly factor encourages the use of renewable energy and reduces the traditional energy used in the optimal mix of energy sources. The proposed model is tested based on demand of energy inspired from real data. These numerical experiments are analyzed to illustrate the impact of encouraging the use of renewable energy sources by introducing the ecofriendly factor and the influence of penalty price and uncertainty in the demand of energy. The results show that in the presence of low penalty price or low uncertainty a large amount of ecofriendly support is needed for using more renewable energy sources in the optimal contract capacity combination.


Motivation
Due to the increase in pollution, the price of oil, and the need of independence from oil, the societies and governments are demanding the reliance on renewable energy sources instead of traditional energy sources. Renewable energy sources are everlasting and nature friendly, while traditional energy sources cause pollution, have huge greenhouse effect on the planet, and are short in supply. The European Union (EU) monitors energy consumption of its countries systematically by defining its energy targets for 2020, 2030 and 2050. To maintain electric stability, availability and different aspects the existing traditional energy such as diesel generators can't be removed, so it is necessary to combine multisources of energy composed of renewable energy sources and nonrenewable energy sources at the same time. In this paper multisources of energy contract capacities are optimized in an efficient and compatible manner.
Practitioners and researchers work on different challenges offered by multisources of energy, such as: location optimization, dimensioning the implementations, optimizing the efficiency with effectiveness taking into account the uncertainties related to the energy generation and the demand. Energy contracting strategies for multisources of energy is one of the important contributions to optimize the costs.

Literature Review
Hybrid energy system composed of wind turbines, solar photovoltaic panels, batteries and auxiliary fuel generator have been mathematically modelled by Gupta et al. [1] and Billionet et al. [2]. Similar problems are also studied with considering different kinds of generators as mentioned by Gupta et al. [1] such as : micro-hydro generator, biomass generator, and biogas generator. The aim of their model is to identify the most economic and appropriate power supply and operating the diesel engine generator under constant output with high efficiency to reduce pollution. The authors considered remote rural areas, and formulated the problem as a mixed integer linear programming with the objective of minimizing the energy costs.
Uncertainties in energy generation and demand had been handled by Prabhu et al. [3], Zhou et al. [4], and Billionet et al. [2]. A two-stage robust approach with constraint generation algorithm is used by Billionet et al. [2], where each sub-problem is reformulated by a mixed-integer, linear program and hence solved by a standard solver. The authors use a polynomial time dynamic programming algorithm for the recourse problem and showed that, in some cases, this algorithm is much more efficient than mixed-integer linear programming approach. Stochastic/probabilistic models for the same problem were considered by Prabhu et al. [3]. To solve this problem Zhou et al. [4] used a two-stage decomposition based solution strategy with genetic algorithm performing the search on the first stage variables and a Monte Carlo method dealing with uncertainties in the second stage.
The problem of energy efficiency versus effectiveness was solved by Prabhu et al. [3]. They optimized energy efficiency considering manufacturing system effectiveness as a constraint, or they optimize manufacturing system effectiveness for a given profile. They used the total available power as a constraint, or they defined a simultaneous optimization of both to get a balanced solution.
The economic efficiency and energy intensity consumption had been analyzed by Bojnec and Papler [5] as determinants of sustainable economic development for 33 selected European countries. The methods used are multivariate factor analyses, regression, and correlation.
The trade off between the costs of renewable energy technologies and the reliability of a multi-sources generation system had been studied by Bilal et al. [6]. They proposed a fast and elitist multi-objective non-dominated sorting genetic algorithm NSGA-II to optimize the use of renewable energy technologies taking into account the annualized cost and the power system reliability in terms of load supplied and including renewable sources in power generation.
Buildings with multisources of energy need optimization to supply its load with less cost and more efficiency. A multisource system technology's proper sizing and energy demands are optimized by Barbieri et al. [7], they used genetic algorithm for minimum energy consumption and net present value, the method is applied on a thirteen-floor tower. To reduce the energy demand, cost and emissions Galvão et al. [8] developed an energy model based on a mixed system of renewable energy. This energy model is a new sustainable standard about energy consumption efficiency of a small hotel building and a relevant contribution to certify the building. For minimum life-cycle costs of meeting the energy demand (power, heating, cooling) of a commercial building, Safaei et al. [9] integrated cogeneration, solar and conventional energy sources. They optimized the investment planning and operating strategies of the energy systems using general algebraic modeling system. The analysis showed varied operating strategies and output levels among cogeneration technologies and the energy systems coupled with them.
Various problems face multisources of energies in microgrids. Achieving better economic and environmental benefits of microgrids under multiple uncertainties in renewable energy resources and loads is problem solved by Wang et al. [10]. The energy production scheduling method is based on robust multi-objective optimization with minimax criterion, they used a hierarchical meta-heuristic solution strategy including multi-objective cross entropy algorithm. Distributed intelligent multi-agent technology is applied by Logenthiran et al. [11] to make the power system more reliable, efficient and capable of exploiting and integrating alternative sources of energy. The simulation results of a power system with distributed resources comprising three microgrids and five lumped loads show that the proposed multiagent system allows efficient management of micro-sources with minimum operational cost. To manage the problematics for different levels of information sharing in a smart grid, Caron and Kesidis [12] proposed a dynamic pricing scheme incentivizing consumers to achieve an aggregate load profile suitable for utilities, and study how close they can get to an ideal flat profile depending on how much information they share. They provide distributed stochastic strategies that successfully exploit this information to improve the overall load profile. For future smart grid, Mohsenian-Rad et al. [13] considered autonomous and distributed demand-side energy management system among users that take advantage of a two-way digital communication infrastructure envisioned. The procedure utilized is game theory and formulation of energy consumption scheduling game. To achieve a stable operation after islanding with minimum fuel cost during the grid-connected opteration, Nam et al. [14] developped economic dispatch problem and the constraints considering reservation for variation in load demand, power outputs of non-dispatchable distributed generators, and stable islanded operation, with flow limits between two different areas. For optimally designing distributed energy resources within cooling and heating power based microgrids, Zhang et al. [15] coupled environmental and economic sustainability in a multi-objective optimisation model using genetic algorithm which integrates the results of a life cycle assessment.
Researchers and practitioners had different attempts for managing energy pricing. A quantitative analysis on real time pricing and the energy market had been studied by Poletti and Wright (2017) [16] for New Zealand introducing different assumptions on the retail market and the shape of the demand function. For scheduling electrical appliances for an individual household, taking into account a grid connected system with a battery and an in-house renewable energy generator, Mitra and Dutta [17] offered dynamic prices as function of the planned consumption and forecasted grid load in the optimization model for minimum electricity bill. To reduce the consumption of electricity in peak periods for industries, Safdarian et al. [18] considered time varying prices by establishing a stochastic model for the medium-term decision making problem faced by a distribution company. The model is formulated as a mixed integer linear programming and using a price elasticity matrix the demand response to time of use prices is captured. To bridge the speed of response gap between suppliers and consumers yet adhere to the principle of marginal cost pricing of electricity, Çelebi and Fuller [19] examined the complementarity programming models of equilibrium. They developped a computable equilibrium model to estimate ex ante time of use prices for retail electricity market. To study the biomass supply contract pricing and policy making in the biofuel industry, Huang and Hu [20] proposed an agent-based simulation model formulated as a mixed integer optimization model. In this model, the agents include a biofuel producer and farmers. Farmers' decision-making is assumed to be profit driven. Due to the dissimilarity between the peak and off-peak hours prices, Sulaima et al. [21] introduced Demand Side Management through Demand Response technique for the modification of the demand profile by implementing different strategies of measures. The objective of this study is to optimize the energy profile for commercial sector, as well as analyze the significance of electricity cost reduction by using Evolutionary Algorithm Meta-heuristic optimization technique in Malaysia. In electricity markets under uncertain price Luo et al. [22] set up a new self-scheduling model based on robust optimization methodology. By using optimal dual theory, the proposed model is reformulated to an ordinary quadratic and quadratic cone programming problems in the cases of box and ellipsoidal uncertainty respectively. The model is tested on IEEE 30-bus system is used to test the new model.
There are attempts for optimizing the total energy costs for machines. For example, two new mathematical models to reduce total energy consumption cost of a single machine manufacturing system are presented by Aghelinejad et al. [23]. The first model improved the formulation of Shrouf et al. [24] problem considering a predetermined jobs sequence. Meanwhile, the second model studies production scheduling on job levels and machine, proposing an optimal sequence for them by minimizing the occurrence number of each machines state, the optimal allocation of these states during the periods minimizes energy costs and jobs within the processing state. The optimal solution of the problem is provided by calculating the shortest path between the first node and last node representing respectively the first and the last periods. Machines consume different amount of energy in function of its state: processing, idle or off state. So for single machine scheduling problems, Aghelinejad et al. [25] proposed a dynamic programming approach to solve these problems by using a finite graph. A non-preemptive single-machine manufacturing environment had been investigated by Aghelinejad et al. [26] to reduce the total costs. They Improved the mathematical formulation of scheduling problem in a predetermined order at machine level to process the jobs. Second, the authors generalized the model to deal simultaneously with the production scheduling at machine level as well as job level. A heuristic algorithm and a genetic algorithm were applied on the second model because the problem is NP hard to provide good solutions in reasonable computational time in an accurate and efficient manner. To minimize the total energy consumption costs over the planning horizon of several energy-oriented single-machine, Aghelinejad et al. [27] studied two cases: constant energy price and increasing energy prices during all the time-slots. Moreover, two versions are investigated: with and without the fixed sequence for the jobs. A multi-item capacitated lot-sizing and scheduling problem had been discussed by Masmoudi et al. [28] in a flow-shop system with energy consideration. They formulated a mixed integer linear programming to solve this issue, and because the capacitated lot-sizing problems are NP-hard, a fix-and relax heuristic is put in place. A single-item capacitated lot-sizing problem in a flow-shop system with energy consideration is also addressed by Masmoudi et al. [29]. In addition to fix-and-relax heuristic method, a genetic algorithm is formed for better quality solutions and to deal with the complexity of the NP-hard problem. In a two stage production system, where a product is manufactured on a machine and delivered to the subsequent production stage in batch shipments, Zanoni et al. [30] proposed an analytical model of this system to minimize the total cost of the production, storing, and energy. To optimally plan saving for energy aware scheduling of manufacturing processes, Bruzzone et al. [31] proposed a mixed integer programming model where the reference schedule is modified to account for energy consumption without modifying the jobs' sequencing and assignment provided by the reference schedule. For sustainable consideration in manufacturing scheduling, Mansouri et al. [32] incorporated energy consumption as an explicit criterion in shop floor scheduling. They explore the potential for energy saving in manufacturing. They analyzed the trade-off between total energy consumption, minimizing makespan, and a measure of service level. To find the Pareto frontier comprised of makespan and total energy consumption, they developed a mixed integer linear multi-objective optimization model. For Sustainable scheduling Gahm et al. [33] developed a research framework for energy-efficient scheduling.
There are different electricity tariffs for residential, commercial and industrial customers. The service types applicable to industrial customers are further classified into low-tension and high-tension. The rate schedules available for high-tension service are based on time of use and maximum demand. This paper focuses on the electricity contract decisions of energy consumers such as high-tension industrial customers.
Most industrial and commercial customers are required to pay for their peak demand, besides the energy they consume. Electric utilities charge them for the highest average demand measured in any 15 or 30 min during their billing period. Billing demand is based on consumers' measured maximal demand and their contract demand with utility companies as per the supply agreement [34].
In mainland China, utilities allow large customers to adjust their contract demand monthly [35], if the peak demand does not exceed the contract demand, a fixed demand charge is levied; on the other hand, if the peak demand exceeds the contract demand, a penalty charge twice as the basic rate is levied. In Southern Africa, the contract demand is determined by the notified maximum demand [35]. Customers can temporarily or permanently increase/decrease their notified maximum demand, and their demand charge is based on the maximum of the measured demand and notified maximum demand. Jemena Electricity Networks Ltd., of Victoria, Australia, also has a similar contract demand reset policy, which allows their customers to permanently/temporarily increase and permanently decrease their contract demand through request. Based on Shanghai, the timely adjustment and notification of the contract demand is a valuable information source for utilities' load forecast and maintenance planning [36]. Along with the development of advanced metering infrastructures and data analysis platform, it is expected that a wider range of customers, to provide more accurate and complete information for smart grid operation and smart city functionalities, can adjust the contract demand more frequently and conveniently.
When it comes to the industrial customer side, matching system requirements with the offers in the energy market is one of the important decisions that must be made [37]. Many industrial customers opt to sign a maximum contracted demand. Such an electricity bill consists of an energy charge and a capacity charge. The energy charge is based on kilowatt-hours, while the capacity charge is based on maximum demand consumed during each time of use period. If the peak demand does not exceed the contract capacity, a fixed capacity charge is levied. On the other hand, if the peak demand exceeds the contract capacity, a penalty charge from two to three times the basic rate is levied [35]. Hence, choosing an excessively low contract capacity will impose high capacity charges, while choosing an excessively high contract capacity may result in an unnecessary basic capacity charge. Therefore, optimal contract capacity decisions have received significant attention from customers with high electricity usage.
In some countries, electricity bills are composed of five primary components: capacity cost, demand of energy cost, power factor adjustment, penalty cost, and expanding construction cost. The energy costs include the active and reactive energy charges but we left the demand of energy out of the context, since it does not affect the optimization problem. Improper contract capacities capacity scheduling would also cause high expanding construction cost when users modify their contract capacities. If customers modify their contract capacities, electricity suppliers would adopt expanding line construction cost schedules. Capacity cost is based on kilowatt hours, with the unit price varying by peak, medium and off-peak and capacity charge is determined by kilowatts per month based on maximum demand (in 15 minutes average) during the time of use period [38]. In our case, we consider only peak contract capacity for simplicity.
There are different contributions to find the optimal contract capacity for energy consumers. To determine the electricity contract capacity for industrial customers in Taiwan Chen et al. [35] used a linear programming approach. To calculate the optimal capacity for peak, semi-peak, and off-peak capacities respectively, Tsay et al. [39] used meta-heuristic evolutionary programming method. The method is applied on Drow-Ing refinery contract, Ho-Jin region contract, Ling-Jan refinery contract, and Da-Liau station contract. To solve the same problem for selection of optimal capacities, for peak, semi-peak, and off-peak capacities Lee et al. [38] used another meta heuristic method called iteration particle swarm optimization. A new index, called iteration best is incorporated into particle swarm optimization to improve solution quality and computation efficiency. Demand contract decision for the taiwanese industries are optimized by Hwang et al. [40]. The techniques employed are cat swarm optimization and particle swarm optimization. This research aims at exploring the benefit on load management options and to provide decision-makers and leaders with useful operation and management strategies as reference. For the assignment of the optimal production planning and energy contract Rodoplu et al. [37] used linear programming in CPLEX software. The optimization model minimize energy and production costs with respect to constraints of energy supplier contract conditions and production systems. As an application, they chose different contract capacities of multisource of energy in three machines for minimum costs and optimal production planning. Optimization of contract capacity setting for industrial consumer with self-owned generating units is a highly discrete and a nonlinear model. To resolve this issue Yang and Peng [41] considered the peak, semi-peak, and off-peak contract capacities in their model. To clear up this muddle the authors used an improved Taguchi method by combining it with particle swarm optimization algorithm. Their paper uses data derived from the SCADA system of a large optoelectronics factory with self-owned generating units in September 2004 to March 2005. When several rates are available in the market Ferdavani et al. [42] proposed new procedure to unfold the contract capacity optimization problem. They include peak contract capacity in their optimization model. The proposed technique gave a better concept for optimizing contract capacity with errors in the forecasted prices or forecasted maximum demand and is faster than the linear programming. In case of variable contract capacity price and uncontracted demand price, using Newton-Raphson method the solution is found within maximum two iterations. To highlight the effectiveness of this method, the proposed approach is performed on the data of various scenarios of a large real electrical user in Singapore. Optimal demand contracting strategy under uncertainty is a complex problem unraveled by Feng et al. [43]. They proved the convexity of the objective function under uncertainty and consequently that there exists one global optimal solution. To calculate the optimal contract value they adopted Newton-Rapson-based numerical method. In the absence of the wide adoption of the energy performance contracting, due to the reasons of its investment-assessing deficiency, high uncertainty, and the complexity profit allocation, Guo et al. [44] applied real option analysis for finding the new metrics of the optimal scheme and incorporated contractual flexibility. The method real option analysis utilizing the binomial tree pricing model evaluates the investment value. It single out the feasible and attractive optimal scheme for the service company and the owner.

Contribution
This paper focuses on the electricity contract capacity considering multisources of energy for electric consumers such as industries, home residence and regions. Industries pay for different electric services, like connection, energy, demand, and reactive power consumption depending on how they use electricity as mentioned by Wang et al. [34]. Most industrial and commercial customers, besides the energy that they consume, are required to pay for their peak demand, and the consumers' contract capacity (CC) and the penalty of the excess peak demand over the contract capacity.
Energy suppliers need to know the capacity demand to plan the generation and the transmission of the energy for better service to the customers. Choosing an excessively high contract capacity may result in an unnecessary basic capacity charge, while choosing an excessively low contract capacity will impose high penalty charges.
For example, when the demand exceeds the upper or lower limits of the contract option, the double power price is charged on the excess quantity as introduced by Rodoplu et al. [37]. In Taiwan, an excess within 10 % of the contract demand is penalized at twice the rate of the contract demand, whereas the excess over 10 % of the contract demand is charged at thrice the rate, and the contract demand can be changed by high-tension industrial customers each month as modeled by Chen and Liao [35].
Multisources of energy contract capacity optimization adding an ecofriendly encouragement factor is introduced by Hamze et al. [45]. In this work, discrete contract capacity of different types are chosen optimally taking into consideration the penalty price for the excess demand over the total combination of contract capacities. The model is solved using linear programming and the effect of the encouragement factor is examined on the total costs (TC) and the total renewable energy contract capacities used. This paper is an extension of our previous work, stochastic features are introduced to the demand of energy considering continuous contract capacities of multisources of energy.
Different types of energy contracts can be proposed by energy suppliers: traditional and renewable energy contracts (solar, wind, biomass ...) as done by Direct Energy electric company in France [46]. The government encourages the use of renewable energy sources by adding bonuses for the renewable energy contracts because of the pollution and for social reasons, it also discourages the use of traditional energy through taxing policy. It is necessary for industries to find the optimal combination of contracts to increase the percentage of green energy used for marketing purposes at the first hand, and to satisfy their peak demand with minimum costs. Moreover, industries them-self support the use of renewable energy contracts. So in general, this support of renewable energy contracts and discouragement of the use of traditional energy contracts is described in this paper by P eco . The objective is to form an optimal mix between the renewable energy energy contract capacity and the traditional non-renewable contract capacity. The more is the percentage of the renewable energy sources used, the more the contracts are considered ecofriendly.
The demand of energy changes from time to time depending on the user and many other factors, so the demand of energy is subjected to uncertainty. This paper addresses a nonlinear model to optimize the energy capacity contracting strategy under demand uncertainty while considering traditional and renewable energy sources as a new concept. This model identifies optimum combination of different types of contract capacities taking into account different factors such as the maximum demand, penalties of over consumption and the ecofriendly price for encouraging the use of renewable energies. The influence of the different values of the ecofriendly price, the penalty prices, and uncertainty on the optimal solution are studied in the presence of the different types of contract capacities. Since the penalty price is previously defined by the producer and the uncertainty is uncontrollable, this approach gives an idea about how much the ecofriendly support should be in the presence of different penalty prices and uncertainties.
The remainder of the paper is organized as follows. In Section 2, the problem statement and the mathematical formulation are described. Section 3 explains the interior point algorithm used in the optimization approach. In Section 4, the model is tested several numerical examples inspired from based on real data. Section 5, summarizes the contribution of this work and introduces some perspectives.

Problem Statement
In this paper, the penalty for excess peak power demand over the total contract capacity values is studied, this excess in power is multiplied by the penalty price P p in ($/ unit of power) . Stochastic features are introduced to the demand, since the demand of energy for the consumer is a future expectation. The types of cost considered are contract capacity cost, penalty cost, and ecofriendly cost. Governments support the use of renewable energy for the consumers, and the consumers such as industries themselves support the use of renewable energy to have green products for marketing purposes. This support is represented by the ecofriendly cost P eco , it gives discount when using more renewable energy, and pays more when using additional traditional energy. Continuous values of the contract capacities are considered, there are traditional contract capacity and different types of renewable energy contract capacities such as solar, thermal, wind, hydroelectric, biomass, bio-fuels, etc. K is the total number of these renewable energy contract capacity types. In the following, a nonlinear model is proposed for this problem. The indexes, parameters and variables of the proposed model are reported in the table of notations.

Mathematical Model
The objective of this model is to find the optimal solution that encourages the use of renewable energy sources. The objective function represents the total cost, it is composed of three parts. The first part is composed of the total contract capacities cost, the second is penalty cost, and the third is the ecofriendly encouragement cost. In this optimization model, the optimal solution of the total contract capacities of the different types should be found to satisfy the power demand of the consumer with minimum cost.
Pren k Ren t,k ) Indexes t(t = 1, ..., T) Index of periods k(l = 1, ..., K) Index of type of renewable energies Variables X t Excess demand at period t Trad t Power capacity of traditional contract at period t Ren k,t Power capacity of renewable contract at period t of type k Parameters T Total number of periods L Total number of contracts K Total number of renewable power contract capacity types P p Penalty price in $/unit of power P eco Ecofriendly price in $/unit of power Ptrad Price of traditional contract in $/unit of power Pren k Price of renewable energy contract of type k in $/unit of power D t Random peak power demand at period t x Possible value ofD t fD t (x) Probability density function ofD t The first part of the objective function shown in Equation (1) represents the sum of the prices of the different chosen contracts. The variables Trad t and Ren t,k will take the value of the optimal combination of traditional and renewable energy contract capacities to be charged with their respective prices Ptrad and Pren k for each period t.
Equation (2) of the second part of the objective function represents the ecofriendly encouragement factor, the more the traditional energy used, the more the consumer pays. On the other hand, the more the consumer uses renewable energy sources, the more there is a discount for the consumer.
The last part of the objective function in Equation (3) represents the penalty cost of the total excess demand estimation for all periods. The penalty price P p in ($/unit of power) is multiplied by the estimated demand exceeding the total optimal contract capacity combination. The probability of excess demand exceeding the total contract capacities are illustrated in Figure 1 considering normal distribution. When the demand is less than total contract capacity combination no penalty is charged. The integration starts from the chosen total contract capacities until the maximum value attainable by the demand of energy. The excess demand estimation is illustrated in Figure 2 as X t .  The objective is to find the optimal combination of the different types of contract capacities of the multisources of energy with minimum costs that satisfy the given demand. The mathematical model is described as following: Subject to: Ren k,t ≥ Ren min ∀ t = 1, ..., T, k = 1, ..., K, Trad t ≤ Trad max ∀ t = 1, ..., T Ren k,t ≤ Ren max ∀ t = 1, ..., T, k = 1, ..., K, Constraints from (5) to (8) are to define the minimum and the maximum value of the traditional contract capacity and the different types of renewable energy contract capacities.

Optimization Approach
In this paper, interior point method is described and analyzed for solving nonlinear constrained optimization problem. The overview of the optimization approach is shown as following: Interior point method Consider the following model: where f : R n −→ R, c : R n −→ R m are smooth functions having derivatives of all orders everywhere in its domain, n is the total number of variables and m is the total number of constraints. The interior point strategy associates logarithmic barrier to the objective function and constraints to solve it and obtain the solution: The following nonlinear system is caused by Karush-Kuhn-Tucker (KKT) conditions of the barrier problem (10): λ ∈ R m and 0 ≤ z ∈ R n represent the Lagrangian multipliers and X = diag(x1, x2, ..., xn). Multiplying the second row of the matrix (11) by X , the system obtained is: This may be viewed as a perturbed KKT system for the original problem (9). The optimality error for the barrier problem is defined based on (12) as: is used to measure the optimality error for the original problem (9).
Primal-dual interior point methods as mentioned by Qiu and Chen (2018) [47] apply Newton's method to the perturbed KKT system and modify step-size so that the inequality (x, z) ≥ 0 is satisfied strictly. A primal-dual linear system is given as: where H is the Hessian of the Lagrangian function and Z = diag(z1, ..., zn). Eliminating d z by d z = −z + µX −1 e − X −1 Zd x , and defining λ + = λ + d λ , we have the iteration It is easy to see that the step generated by this system coincides with the solution of the following primal-dual quadratic programming subproblem Step computation of the Algorithm 1 is based on this model, λ + = λ + d λ , x + = x + d x , and z + = z + d z . The parameter K needs to be selected and the barrier parameter µ should be updated according to the procedure in literature [48][49][50] for fast convergence. So the Algorithm 1 is expressed as the following:

Algorithm 1: Interior Point Algorithm initialization:
Choose the starting point (x 0 , λ 0 , z 0 ) ; Select barrier parameter µ 0 > 0; parameter end Algorithm 1 can solve linear and nonlinear convex optimization problems. In this study the objective function is the total cost and the constraints are the boundaries of the different types of contracts. To prove the convexity of the problem we should study the hessian matrix of the total cost since it is a multivariable function. The total cost is a function of Traditional and Renewable energy contract capacity variables. The determinant of the hessian matrix and the diagonal components of the matrix should be positive. Based on Leibniz integral rule for partial derivation of the total cost ∇ 2 TC is calculated and the followings are obtained: Ren k )) (16) Ren k ) (20) The hessian matrix becomes as following: The components of the hessian matrix are all the same so the determinant is equal to 0, and the diagonal components (17) and (19) are strictly positive since P p > 0 and fD(Trad + ∑ K k=1 Ren k )) > 0 so the problem is convex.

Experiments Design
The model and numerical experiments design are inspired from the work of Feng et al. [43]. In their work, the authors tried to find the optimal contract under uncertainty in case of one period demand and one kind of energy source for different probability distributions. The proposed model is tested based on adopted real data representing the monthly energy consumption of Grand-Est Region (in France) [51] for year 2018. The data may follow any type of distribution, here the data are assumed to follow a normal distribution density function: The normal distribution density function is used because its parameters directly represent the mean and the variance, by changing the variance the effect of uncertainty will be tested directly. In this study the contracts are assumed to be constant in all periods. So the following constraints are added: Ren k,t+1 = Ren k,t ∀ t = 1, ..., T − 1, k = 1, ..., K, In this case, two renewable energy contract capacities are considered (solar and wind energies K = 2). For the set of data the problem is optimized to find the best combination of different types of contract capacities. The effect of the penalty price on the optimal solution in one type of traditional energy is well known, but in the presence of multi types of renewable energy contracts, it is needed to be discovered. To study the influence of the ecofriendly price, penalty price and uncertainty it is necessary to change the values of P eco , P p , and σ as shown in the graph of Figures 3-11. The values of the parameters of the model are assumed as shown in Table 1. These parameters are presented by the prices in $/MW, the boundaries are the maximum and minimum values in MW achievable of the different types of contract capacities of traditional, solar and wind energy types. When studying the effect of ecofriendly price P eco and penalty price P p a certainty case of σ = 0 is considered and the values of P p are shown and described in Table 2  Price value in ($/MW)  7500  8070  8750 18,000 30,000 In comparison with the contract capacity price P p1 < P Trad P Trad < P p2 < P Sol P Sol < P p3 < P Wind P p4 > P Wind P p5 > P p4 Calculation with respect to contract capacity price - PSol +PWind 2 2 * P Wind 3 * P Wind Meanwhile when analyzing the effect of ecofriendly price P eco and uncertainty, the value of the penalty price remains constant of P p = 18,000 $/MW equal two times of the highest price of contract capacities P Wind . The demand of energy are stochastic and they follow a normal distribution of average value presented in the Table 3, the standard deviation is changed to study the effect of uncertainty on the optimal solution. The values of Table 3 are the energy demand in MW of the Grand-Est region of France for the year 2018. The value of the standard deviation of the data in Table 3 is calculated to be σ = σ 3 = 596.9702, so the values of the standard deviations representing the uncertainty studied are multiple of σ 3 as shown in Table 4, the different values of σ are tested as shown in Figures 8-11. Table 3. Energy demand of Grand-Est for 2018 [51].

Discussion
The problem is solved on MATLAB software using interior point Algorithm 1. The initial point for Algorithm 1 used is equal to the lower bound Trad min = 500 MW, Sol min = 250 MW, and Wind min = 250 MW. The comparison is done between the rate of change of the different types of contract capacities. The effect of P eco and P p on the optimal solution are studied by changing their values.
Figures 3-7 present the variation of the optimal contract capacity types, total cost, and total excess demand with respect to P eco and P p in case of certainty σ = 0, the effect of P p on the optimal solution is examined. For the values of penalty price less than the traditional energy price which is the minimum between the different types of contract capacities, the optimal solution is the minimum value of the contract capacities Trad min = 500 MW, Sol min = 250 MW, and Wind min = 250 MW. That is because the penalty is cheap so even if the demand exceeds the contract capacities too much the price of the penalty would be cheaper than that of the contract capacities.
As the P p increases above the minimum price of the contract capacities, in this case the traditional energy, the optimal solution of the contract capacities increases starting from the traditional energy from 500 MW until reaching 3000 MW as shown in Figure 3, but when the P eco increases the optimal traditional energy drops to minimum 500 MW, since P eco discourages the use of traditional energy by increasing the cost.
As the P p increase the optimal solar contract capacity augments from 250 MW to 1120 MW for P eco = 0. As the P eco rises the optimal solar contract capacity rises from 250 MW till reaching 2200 MW in one step since it is cheaper than the wind contract price as presented in Figure 4. The wind energy contract capacity increases continuously with P eco since it is the most expensive contract capacity as can be seen in Figure 5. As the P p augments the optimal wind energy contract capacity augments more rapidly with P eco till reaching saturation Wind max = 2200 MW, to compensate the high value of the penalty price subjected on the excess demand.
When the P eco increases the total cost decreases in general as shown in Figure 6. As the P p increases the total cost increases and have the same sense of variation with P eco , but for P eco ≥ 6000 $/MW the total cost is approximately the same for all values of P p , that is because the optimal renewable energy contracts become higher for their cheap price and cover all the expected demand, so there will be no excess demand to pay penalty over and the total cost becomes approximately the same.
The total excess demand decreases as the P eco increases for all values of P p as presented in Figure 7. The P eco increases the optimal renewable contract capacities that cover the estimated demand. As the P p rises the total excess demand reduces, the optimal contract capacity combination ameliorate to minimize the excess demand so the user isn't charged by a high penalty. Figures 8 and 9 represent the variation of the different types of contract capacities with respect to P eco and σ considering constant penalty price of P p = 18,000 $/MW equal two times that of the wind contract capacity price. For the study of the uncertainty's influence, when the P eco rises the optimal traditional energy contract drops from maximum 3000 MW to minimum 500 MW, while the first to increase is the optimal solar energy contract from 250 MW to 2200 MW since it is cheaper. The values and the sense of variation of the traditional and solar contract capacities with respect to P eco are the same for all values of σ as viewed in Figure 8.
In Figure 9, as the uncertainty increases the increase of wind energy with the P eco reaches saturation of Wind max =2200 MW much faster, since the increase of uncertainty gives higher value of the expected demand, so higher values of contract capacities are needed to compensate the high expected demand.
Based on the Figure 10, the total costs decrease obviously with the P eco since it is encouraging the use of the renewable energy after all, but as the value of σ rises the total costs increase, that is due to the high value of the expected excess demand, so the user is charged by the high contract capacity to compensate for the excess demand or the penalty for the expected demand exceeding the total contract capacities.        The total excess demand (TED) is calculated by the difference of the average demand exceeding the total capacity as following: Ren k ). The total excess demand decreases with P eco since there is enough contract capacity when encouraging the use of renewable energy covering the expected demand. With the varaition of uncertainty the total excess demand decreases more rapidly for the reason of the increase of the expected demand with uncertainty, so the total optimal contract capacities increases so that the user isn't charged by high penalty as seen in Figure 11. Figure 11. Variation of the total excess Demand with respect to P eco and σ (own results).
The Table 5 shows the variation of total contract capacities, total cost, and total excess demand with respect to P eco , P p , and σ, the positive sign signify increase, and the negative sign signify decrease. Penalty price is subjected to the government laws, and the uncertainty depends on the user way of consumption, so these results show that in case of high penalty low values P eco are needed to increase the use of renewable energy. When there is enough data about the demand of energy high amount of support of P eco is needed to use renewable energy because the expected demand will be low. Table 5. Variation of contract capacities, total cost, and total excess demand with respect to P eco , P p , and σ (own results).

Different Distributions
The problem is tested on two other different probability density functions, gamma and log-normal distributions as described in appendices A and B, to illustrate the generality of the model and compare the results between them. The parameters of the gamma and log-normal distributions are calibrated to have the same mean and standard deviations tested. The results of the optimal wind contract capacity and total cost are presented in Figures 12 and 13 respectively in case of gamma distribution, for the case of the log-normal distribution they are presented in Figures 14 and 15. The total excess demand in case of gamma distribution and log-normal distribution are both shown in Figure 16.     The optimal traditional and solar wind contract capacities in case of gamma and log-normal distributions are approximately the same as in Figure 8. The results in general obtained from the gamma and log-normal distribution are approximately the same as the results obtained from the normal distribution. But in case of total excess demand for σ it is the same for the different distributions, but for 2σ and 3σ the log-normal has the highest total excess demand, second is the gamma distribution, then the normal distribution as presented in Figure 16. So the influence of uncertainty on the expected excess demand changes from one type of probability distribution to other depending on its form, but it is not significant enough to change the the optimal contract capacities and the total cost.

Conclusions
In this paper we have formulated a model to optimize the total cost for choosing the best contract capacity combination of multisources of energy. To solve the model we used interior point algorithm taking into account the cost of the different types of contract capacities, the penalty price of the excess demand, ecofriendly encouragement for the use of renewable energy, and uncertainty in the energy demand. The multisources of energy contract capacities are composed mainly of traditional and one or more types of renewable energy such as solar and wind in the case studied.
The penalty price increases the optimal contract capacity starting from the contract with minimum price, until the penalty price exceeds the contract with the maximum price, the contract capacities combination in total increase, that is for the demand not to exceed the total contract capacities and be paid by large penalty costs. But in the presence of the ecofriendly factor as it increases with the increase of the penalty price, the renewable energy contract capacities increases and traditional energy contract capacity decreases in the combination of contract capacity, this combination increases in total with the penalty price. So in the situations having high penalty prices it is not necessary to have large ecofriendly support to encourage the use of renewable energy in the contract capacity. But for low penalty price a great amount payment is used to have more renewable energy in the contract capacity combination.
Uncertainty is studied for the demand of energy and changing the ecofriendly factor at the same time, the stochastic features are examined in different probability distributions such as normal, gamma, and log-normal probability distributions. As the uncertainty increases the expected excess demand gets higher and eventually the contract capacities of all types including renewable energy sources and total cost rise while the total excess demand reduce. Therefore, when having lack of data concerning the demand of energy or high uncertainty for any reason, a small ecofriendly assistant is enough to have more renewable energy in the contract capacity combination but with high total cost and vice versa. The total cost increases with the penalty price and uncertainty which are uncontrollable parameters, but P eco is controllable and the total cost decreases with respect to P eco , so this study shows in the presense of low or high values of penalty price and uncertainty what value of P eco to choose. Great thanks for the three anonymous reviewers for there notes to improve the quality of this approach.
In this paper, we proposed an optimized energy contracting strategy tacking into account the uncertainties related to the energy demand. We assumed these uncertainties represented by any general probability distribution, from the optimization modeling perspective; we can improve current work by investigating robust pricing policies under demand uncertainty, which guarantee a certain level of performance for possible scenarios, even for the worst-case scenario. Also, uncertainty in other parameters can be another extension. For example, a direct extension of this work is to consider the uncertainty related to the availability of renewable energies. It could be interesting to consider both demand and production uncertainties to find the set of contracts of the different types of renewable energy for the consumer to choose.
Author Contributions: All authors have worked on this manuscript together and all authors have read and approved the final manuscript.

Acknowledgments:
The authors are grateful to UTT in France, FEDER, and the Doctoral School of Science and Technology in the Lebanese University for their financial support. The authors would like to thank the three anonymous reviewers for their valuable and helpful remarks.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:  The parameters a and b of the gamma distribution function have been calibrated in the manner to have the same mean µ of the demand presented in the Table 3 and standard deviations studied in Table 4. The parameters a and b of the log-normal distribution function have been calibrated in the manner to have the same mean µ of the demand presented in the Table 3 and standard deviations studied in Table 4.