A Bilevel Approach to Optimal Price-Setting of Time-and-Level-of-Use Tariffs
Mathieu Besançon, Miguel F. Anjos, Luce Brotcorne, Juan A. Gomez-Herrera
EEUROGEN 2019 September 12-14, 2019, Guimarães, Portugal
A Bilevel Approach to Optimal Price-Setting of Time-and-Level-of-Use Tariffs
Mathieu Besançon
Department of Mathematics & Industrial Eng., Polytechnique Montréal & INRIA Lille & Centrale LilleINRIA Lille, 40 avenue Halley, 59650 Villeneuve-d’Ascq, FranceEmail: [email protected]
Miguel F. Anjos*
School of Mathematics, University of EdinburghJames Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh EH9 3FD, UKEmail: [email protected]
Luce Brotcorne
INOCS, INRIA LilleINRIA Lille, 40 avenue Halley, 59650 Villeneuve-d’Ascq, FranceEmail: [email protected]
Juan A. Gomez-Herrera
GERAD & Department of Mathematics & Industrial Eng., Polytechnique Montréal2900 Boulevard Edouard-Montpetit, Montréal, QC H3T1J4, CanadaEmail: [email protected]
Summary
Time-and-Level-of-Use (TLOU) is a recently proposed pricing policy for energy, extending Time-of-Use with the additionof a capacity that users can book for a given time frame, reducing their expected energy cost if they respect thisself-determined capacity limit. We introduce a variant of the TLOU defined in the literature, aligned with the supplierinterest to prevent unplanned over-consumption. The optimal price-setting problem of TLOU is defined as a bilevel,bi-objective problem anticipating user choices in the supplier decision. An efficient resolution scheme is developed,based on the specific discrete structure of the lower-level user problem. Computational experiments using consumptiondistributions estimated from historical data illustrate the effectiveness of the proposed framework.
Keywords:
Demand response, electricity pricing, bilevel optimization, Time-and-Level-of-Use
The increasing penetration of wind and solar energy hasmarked the last decades, bringing a decentralization andhigher stochasticity of power generation, yielding newchallenges for the operation of electrical grids. Advances incommunication technologies have enabled both schedulingof smart appliances and seamless data collection andexchange between different entities operating on powernetworks, from generators to end-consumers. Usingsuch capabilities, agents can make decisions based onprobability distributions of different variables, possiblyconditioned on known external factors, such as the weather.As a promising lead to tackle these challenges, Demand
This extended abstract was presented at, and included in the programof, the 13th EUROGEN conference on Evolutionary and DeterministicComputing for Industrial Applications.
Response (DR) has attracted an increasing interest in thepast decades. Instead of compensating ever-increasingfluctuations of renewable production with conventionalpower generation, this approach consists in leveraging theflexibility of some consuming units by reducing or shiftingtheir load. The Time-and-Level-of-Use system presentedin this paper is primarily a price-based Demand Responseprogram based on the classification found in the literature, although the self-determined consumption limit creates anincentive for respecting the commitment on the part of theuser. Its design makes opting out of the program a naturalextension and thus requires less coupling and interactionbetween users and suppliers to work. In a 2012 report,the US Federal Energy Regulatory Commission identifiedthe lack of short-term estimation methods as one of thecritical barriers to the effective implementation of DemandResponse. In a related approach, an incentive-based1 a r X i v : . [ m a t h . O C ] D ec UROGEN 2019 September 12-14, 2019, Guimarães, PortugalDemand Response program is developed with userssignaling both a predicted consumption and a reductioncapacity, with the supplier selecting reductions randomly.Both the program developed by the authors and the TLOUpolicy require a signal sent from the user to the supplier.Bilevel optimization has been used to model and tackleoptimization problems in energy networks. Applicationson power systems and markets are also mentionedin several reviews on bilevel optimization.
Itallows decision-makers to encapsulate utility-generators,user-utility interactions, or define robust formulations ofunit commitment and optimal power flow. Some recentwork focuses on bilevel optimization for demand-responsein the real-time market. Multiple supplier settingsare also considered, for instance in the determination ofTime-of-Use pricing policies, where competing retailerstarget residential users.The Time-and-Level-of-Use (TLOU) pricing scheme was built as an extension of Time-of-Use (TOU),targeting specifically the issue with current large-scaleDemand-Response programs identified in the FERCreport. The optimal planning and operation of a smartbuilding under this scheme was developed, taking theTLOU settings as input of the decision. In this work, weconsider the perspective of the supplier determining theoptimal parameters of TLOU.We integrate the optimal user reaction in the supplierdecision problem, thus modelling the interaction as aStackelberg game solved as a bilevel optimization problem.The specific structure of the lower-level decision isleveraged to reduce the set of possibly optimal solutionsto a finite enumerable set. Through this reduction, thelower-level optimality conditions are translated into a set oflinear constraints. The set of optimal pricing configurationsis derived for different distributions corresponding todifferent time frames. These options can be computed inadvance and then one is chosen by the supplier ahead ofthe consumption time to create an incentive to book andconsume a given capacity. We show the effectiveness ofthe method with discrete probability distributions computedfrom historical consumption data.
The TLOU policy was initially developed from theperspective of smart consumption units which can bea building, an apartment or a micro-grid monitoring itsconsumption and equipped with programmable consumingdevices. It is a pricing model for energy built upon theTime-of-Use (TOU) implemented by several jurisdictionsand for which the energy price is fixed by intervalsthroughout the day. TLOU extends TOU by allowing usersto self-determine and book an energy capacity at each timeframe depending on their planned requirements; and bydoing so, they provide the supplier with information on the intended consumption. We will refer to the capacityas the amount of energy booked by the user for a given timeframe, following the same terminology as the referenceimplementation. Energy prices still depend on the timeframe within the day, but also on the capacity booked bythe user. This pricing scheme is applied in a three-phaseprocess:1. The supplier sends the pricing information to the user.2. The user books a capacity from the supplier for thetime frame before a given deadline. If no capacity hasbeen booked, the Time-of-Use pricing is used.3. After the time frame, the energy cost is computeddepending on the energy consumed x and bookedcapacity c : • If x ≤ c , then the applied price of energy is π L ( c ) and the energy cost is π L ( c ) · x . • If x > c , then the applied price of energy is π H ( c ) and the energy cost is π H ( c ) · x .In the model considered in this paper, only the firsttwo steps involve decisions from the agents, making thedecision process a Stackelberg game given the sequentialityof these decisions. The price structure is composed of threeelements: a booking fee K , a step-wise decreasing function π L ( c ) representing the lower energy price and a step-wiseincreasing function π H ( c ) representing the higher energyprice. The steps of the lower and higher price functions aregiven at different breakpoints: { c L , c L , c L , ... } = C L & { π L , π L , π L , ... } = π L { c H , c H , c H , ... } = C H & { π H , π H , π H , ... } = π H π L ( c ) will refer to the function of the capacity and π Lj tothe value of the lower price at step j . In the initial versionof the pricing, the energy consumed above the capacity ispaid at the higher tariff while the rest is paid at the lowertariff.Considering that most power systems try to preventover-consumption and unplanned excess consumption, weintroduce a variant where the whole energy consumed ispaid at the lower tariff if it is less or equal to the capacity,and at the higher tariff otherwise, see Equation (1). In otherwords, if the consumption over the time frame remainsbelow the booked capacity, the effective energy price isgiven by the lower tariff curve; if the consumption exceedsthe booked capacity, the energy price is given by the highertariff. The total cost for the user associated with a bookedcapacity c and a consumption X t for a time frame t is: C ( c t ; X t ) = (cid:40) K · c t + π L ( c t ) · X t , if X t ≤ c t , K · c t + π H ( c t ) · X t otherwise. (1)2UROGEN 2019 September 12-14, 2019, Guimarães, PortugalThe load distribution X t is a random variable; both userand supplier make decisions on the expected cost over theset of possible consumption levels Ω , given as: C ( c t ) = E Ω [ C ( c t ; x ω )] (2)where E Ω is the expected value over the support Ω of theprobability distribution.Furthermore, π L ( c = ) = π H ( c = ) = π ( t ) , with π ( t ) the Time-of-Use price at the time frame of interest t . This property allows users to opt-out of the programfor some time frames by simply not booking any capacity.TLOU is designed for the day-ahead market, where boththe pricing components and the capacity are chosen aheadof the consumption time. It can however be adjusted toother markets or intra-day settings. The entity definingthe TLOU pricing can also extend the possible settings bydecoupling the time spans for one price setting choice fromthe time frames for capacity booking. For instance, a pricesetting can be chosen by the utility for the week, while thebooking of capacity occurs the day before the consumption.TLOU offers the user the possibility to reduce theircost of energy by load planning, and offers the supplierthe prospect of improved load forecasting, becauseunder-consumption is paid by the excess booking cost whileover-consumption is paid by the difference between higherand lower tariffs.We illustrate this phenomenon with a supplier decision ( K , π L , π H ) on Figures 1 and 2 using the relative cost: x (cid:55)→ C ( c ; x ) x (3)for different values of c . For c >
0, the fixed cost K · c makes low consumption levels expensive per consumedunit, while the transition from lower to higher price makesover-consumption more expensive than the baseline. In theexample, c = . c = . c = . P r i c e $$ / k W h $ | B oo k i n g f ee $ Lower and higher prices Booking fee = 0.3
K cL / H
Figure 1: Example of TLOU pricing R e l a t i v e c o s t $ / k W h Relative cost C(c;x) / x c = 0.0c = 1.5c = 3.0c = 3.5
Figure 2: Relative cost of energy vs consumption fordifferent capacitiesThe supplier first builds their set of options based onprior consumption data. In the proposed method, theprior discrete distribution used can be conditioned on someindependent variables if they are known and influencethe consumption (e.g. forecast external temperature orday of the week). They can then pick a pricing settingfor a given day based on generation-side considerationsand constraints, including the option to stay at a flatTime-of-Use tariff for some or all time frames.
The supplier wishes to determine an optimal set of pricingoptions at any capacity level. In the model developed in thissection, we consider a discrete probability distribution witha finite support, derive some properties of the cost structure3UROGEN 2019 September 12-14, 2019, Guimarães, Portugalwhich we then leverage to formulate the optimizationproblem in a tractable form.At any capacity level, the decision process of thesupplier involves two objectives, the expected revenue fromthe tariff and the guarantee of an upper bound on theconsumption. The expected revenue is given by: C ( c ) = K · c + ∑ ω ∈ Ω − ( c ) x ω p ω π L ( c ) + ∑ ω ∈ Ω + ( c ) x ω p ω π H ( c ) (4)with any capacity booked defining a partition of the set ofscenarios: Ω − ( c ) = { ω ∈ Ω , x ω ≤ c } (5) Ω + ( c ) = { ω ∈ Ω , x ω > c } (6)The function C ( c ) is minimized by the user with respectto their decision c . It is non-linear, non-smooth anddiscontinuous because of the partition of the scenarios by c and the transition between steps of the pricing curves π L ( c ) , π H ( c ) . Both the user and supplier problems are thusintractable with this initial formulation. Proposition 3.1shows that only a discrete finite subset of capacity valuesare candidates to optimality for the user. Proposition 3.1.
The optimal booked capacity for a userat a time frame t belongs to a discrete and finite set ofcapacities S t , defined as:S t = { } ∪ C L ∪ Ω t (7) with Ω the set of consumption scenarios.Proof. The user objective function is the sum of thebooking cost and the expected electricity cost. The bookingcost is linear in the booked capacity, with a positiveslope equal to the booking fee. The expected electricitycost is piecewise constant in the booked capacity, withdiscontinuities at steps of both of the price curves becauseof the π L and π H prices and at possible load levels becauseof the transfer of a load from Ω + to the Ω − set. This canbe highlighted using the indicator functions associated witheach of the two sets: − ( ω , c ) = (cid:40) , if x ω ≤ c ,0 otherwise. (8) + ( ω , c ) = − − ( ω , c ) (9)The expression of the user expected cost becomes: C ( c ) = K · c + ∑ ω ∈ Ω x ω · p ω · (cid:0) π L ( c ) · − ( ω , c ) + π H ( c ) · + ( ω , c ) (cid:1) (10)The sum of the two terms is therefore piecewiselinear with a positive slope. On any interval between thediscontinuity points, the optimal value lies on the lower bound, which can be any point of C L , C H , Ω or 0.Furthermore, let C ( c ) be the user cost for a bookedcapacity and ¯ c such that ¯ c ∈ C H and ¯ c / ∈ { } ∪ C L ∪ Ω . Thehigher tariff levels are monotonically increasing. Let ε > π H ( ¯ c − ε ) = π n , π H ( ¯ c + ε ) = π n + > π n , π L ( ¯ c − ε ) = π L ( ¯ c + ε ) = π Lm , (cid:64) x ω , ω ∈ Ω s.t. ¯ c − ε ≤ x ω ≤ ¯ c + ε . The last condition guarantees there is no load value in the [ ¯ c − ε , ¯ c + ε ] interval and can also be expressed in terms ofthe two load sets split by the capacity: Ω + ( ¯ c − ε ) = Ω + ( ¯ c + ε ) and Ω − ( ¯ c − ε ) = Ω − ( ¯ c + ε ) . Then if such ε exists, we find that: C ( ¯ c − ε ) = K · ( ¯ c − ε ) + ∑ ω ∈ Ω − ( ¯ c ) π Lm · x ω + ∑ x ∈ Ω + ( ¯ c ) π Hn · x ω , C ( ¯ c + ε ) = K · ( ¯ c + ε ) + ∑ ω ∈ Ω − ( ¯ c ) π Lm · x ω + ∑ ω ∈ Ω + ( ¯ c ) π Hn + · x ω , C ( ¯ c + ε ) − C ( ¯ c − ε ) = ε K + ∑ ω ∈ Ω + ( ¯ c ) ( π Hn + − π Hn ) · x ω , C ( ¯ c + ε ) − C ( ¯ c − ε ) > . The discontinuity on any x ∈ C H is therefore alwayspositive and cannot be a candidate for optimality. It followsthat optimality candidates are restricted to the set S = { } ∪ C L ∪ Ω .Proposition 3.1 means we can replace the continuousdecision set of capacities with a discrete set that can beenumerated.The guarantee of an upper bound on the consumption G corresponds to the incentive given to the user againstconsuming above the considered capacity. It is the secondobjective, given by the difference in cost at the capacity,which is the immediate difference in total cost at thetransition from lower to higher tariff: G ( c , π L , π H ) = c t · (cid:0) π H ( c ) − π L ( c ) (cid:1) . (11)The supplier needs to include the user behavior andoptimal reaction in their decision-making process, whichcan be done by a bilevel constraint: c t ∈ arg min c C ( c ) . (12)The user thus books the least-cost option at each timeframe, given the corresponding probability distribution.Given the finite set of optimal candidates S t defined in 3.1,this constraint can be re-written as: C ( c t ) ≤ C ( c ) ∀ c ∈ S t . (13)4UROGEN 2019 September 12-14, 2019, Guimarães, PortugalIf multiple choices of c yield the minimum cost, thechoice of the user is not well-defined. The supplierwould want to ensure the uniqueness of the preferredsolution by making it lower than the expected cost ofany other capacity candidate by a fixed quantity δ > δ ). It is a parameter of thedecision-making process, estimated by the supplier. Thelower-level optimality constraint is then for a preferredcandidate k : C ( c tk , K , π L , π H ) ≤ C ( c tl , K , π L , π H ) − δ ∀ l ∈ S t \ k . (14) In order to obtain regular price steps, the contract betweensupplier and consumer can include further constraints onthe space of pricing parameters. We include three types ofconstraints: lower and upper bounds on the booking fee K ,minimum and maximum increase at each step of the higherprice and minimum and maximum decrease at each stepof the higher price. All these can be expressed as linearconstraints, and we gather them under the constraint set: ( K , π L , π H ) ∈ Φ . (15) The model is defined for each of the capacity candidatesand is thus noted P tk for candidate k and time frame t :max K , π L , π H (cid:0) C ( c tk , K , π L , π H ) , G ( c tk , π L , π H ) (cid:1) (16) ( K , π L , π H ) ∈ Φ (17) C ( c tk , K , π L , π H ) ≤ C ( c lt , K , π L , π H ) − δ ∀ l ∈ S t \ k , (18)where C ( c tk , K , π L , π H ) = K · c tk + ∑ ω ∈ Ω − t ( c tk ) x ω p ωπ Lk + ∑ ω ∈ Ω + t ( c tk ) x ω p ω π Hk , (19) G ( c tk , π L , π H ) = c tk · (cid:0) π H ( c ) − π L ( c ) (cid:1) . (20) The proposed model was implemented and tested usinghistorical consumption data measured on a pilot house.The instantaneous consumption was measured everytwo minutes during 47 months on a residential buildingby the energy supplier and grid operator EDF. Sincethe focus is the energy consumed within a given timeframe, the instantaneous power can be averaged foreach hour, yielding the energy in kW · h and avoidingissues of missing measurements in the dataset. Datapreprocessing, density estimation and discretization, and visualization were performed using Julia, matplotlib and KernelDensity.jl . The construction and optimizationof the model were carried out using CLP from theCOIN-OR project as a linear solver and JuMP. Forall experiments where it is not specified, an inertia of δ = .
05 has been applied. For every time frame and for allcapacity candidates, the bi-objective supplier problem withobjectives ( C , G ) is solved with the ε -constraint methodimplemented in MultiJuMP.jl . In all cases, the objectivesare found to be non-conflicting, in the sense that theutopia point of the multi-objective problem is feasible andreached. This implies that a lexicographic multi-objectiveoptimization solves the problem and reaches the utopiapoint, but does not guarantee that this holds for all problemconfigurations.Figure 3 shows the number of options computed at eachhourly time frame and Figure 4 the capacity level in kW · h of each option. The option of a capacity level of 0 isalways possible. Two examples of TLOU settings obtainedare presented in Figures 5 and 6.The expected cost for the user of booking any capacity,given the price setting provided in Figure 5 is shown Figure7. The most notable result is that in all cases tested,the utopia point, defined as the optimal value of the twoobjectives optimized separately, is reachable. This resultis conceivable given that the supplier decision is takenin a high-level space, allowing multiple solutions to beoptimal with respect to the revenue. In order to ensure theguarantee-maximizing optimal solution, a two-step processlexicographic multi-optimization procedure is used:1. Solve the revenue-maximizing problem to obtain themaximal reachable revenue v .2. Solve the guarantee-maximizing problem, whileconstraining a revenue C ( · ) ≥ v .Figures 8 and 9 present a TLOU configurationoptimized for costs only and for cost and guarantee.All the models solved are linear optimization problemswith a fixed number of variables and a number ofconstraints growing linearly with the number of scenariosconsidered. However, all these constraints are of type: C ( c tk , K , π L , π H ) ≤ C ( c tl , K , π L , π H ) − δ ∀ l ∈ S t \ k . (21)This is equivalent to: C ( c tk , K , π L , π H ) ≤ min l ∈ S t \ k C ( c tl , K , π L , π H ) − δ . (22)Therefore, at most one of the l will be active with anon-zero dual cost. The method can thus be scaled to agreater number of scenarios by adding these constraints onthe fly.With the current discretized distributions containingbetween 5 and 10 scenarios, the mean and median times5UROGEN 2019 September 12-14, 2019, Guimarães, Portugal N u m b e r o f o p t i o n s Number of capacity options at different hours, = 0.05
Figure 3: Number of options at different hours O p t i o n ( k W h ) Capacity options at different hours, = 0.05
Figure 4: Capacity levels of options at different hoursto compute the whole solutions for all candidate capacitiesare below of a second. These metrics are obtained usingthe BenchmarkTools.jl package. A study of the influence of the δ parameter issummarized Figure 10. For any hour, there always exists amaximum value δ max above which it becomes impossible tomake a solution better for the lower-level than the baselinewith a difference greater than δ max . TLOU is designed to price energy across time and toreflect varying costs and requirements from the generationside. Defining two objectives for the supplier, we builtthe set of cost-optimal price settings maximizing theguarantee in a lexicographic fashion. Computations ondistributions built from real data show the effectivenessof the method, requiring a low runtime to compute ( $ / k W . h ) Tariff curves c = 3.47 kW . h K = 0.1949$/kW.h booked capacity
Figure 5: TLOU price settings to incentive for capacity c = . kW · h ( $ / k W . h ) Tariff curves c = 4.51 kW . h K = 0.1531$/kW.h booked capacity
Figure 6: TLOU price settings to incentive for capacity c = . kW · h the set of solutions. Future research will considercontinuous probability distributions of the consumption andthe price-setting problem with multiple users. This work was supported by the NSERC Energy StorageTechnology (NEST) Network.
References [1] Albadi, M. H. and El-Saadany, E. F. A summaryof demand response in electricity markets.
ElectricPower Systems Research (11), 1989–1996 (2008).[2] US Federal Energy Regulatory Commission (FERC).Assessment of Demand Response and AdvancedMetering, , 12 (2012).6UROGEN 2019 September 12-14, 2019, Guimarães, Portugal c o s t ( $ ) User utility c = 3.47 kW . h candidatesoptimumequilibrium limit Figure 7: Expected cost for any capacity booked under theprice settings of Figure 5[3] Vuelvas, J., Ruiz, F., and Gruosso, G. Limitinggaming opportunities on incentive-based demandresponse programs.
Applied Energy , 668 – 681(2018).[4] Dempe, S., Kalashnikov, V., Prez-Valds, G. A.,and Kalashnykova, N.
Bilevel ProgrammingProblems: Theory, Algorithms and Applications toEnergy Networks . Springer Publishing Company,Incorporated, (2015).[5] Colson, B., Marcotte, P., and Savard, G. An overviewof bilevel optimization.
Annals of OperationsResearch (1), 235–256 (2007).[6] Shi, X. and Xia, H. S. Model and interactive algorithmof bi-level multi-objective decision-making withmultiple interconnected decision makers.
Journal ofMulti-Criteria Decision Analysis (1), 27–34 (2001).[7] Fernández-Blanco, R., Arroyo, J. M., andAlguacil, N. On the Solution of Revenue- andNetwork-Constrained Day-Ahead Market ClearingUnder Marginal Pricing Part I: An Exact BilevelProgramming Approach. IEEE Transactions onPower Systems (1), 208–219 jan (2017).[8] Fernández-Blanco, R., Arroyo, J. M., andAlguacil, N. On the Solution of Revenue- andNetwork-Constrained Day-Ahead Market ClearingUnder Marginal Pricing Part II: Case Studies. IEEETransactions on Power Systems (1), 220–227 jan(2017).[9] Street, A., Moreira, A., and Arroyo, J. M. Energyand reserve scheduling under a joint generation andtransmission security criterion: An adjustable robustoptimization approach. IEEE Transactions on PowerSystems (1), 3–14 Jan (2014). ( $ / k W . h ) Tariff curves c_{13} = 3.3kW.h K = 0.3237$/kW.h booked capacity
Figure 8: TLOU price settings for revenue maximizationonly ( $ / k W . h ) Tariff curves c_{13} = 3.3kW.h K = 0.2941$/kW.h booked capacity
Figure 9: TLOU price settings after lexicographicoptimization[10] Bahrami, S., Amini, M. H., Shafie-khah, M., andCatalao, J. P. S. A Decentralized RenewableGeneration Management and Demand Response inPower Distribution Networks.
IEEE Transactions onSustainable Energy (4), 1783–1797 (2018).[11] Yang, J., Zhao, J., Wen, F., and Dong, Z. Y.A framework of customizing electricity retailprices. IEEE Transactions on Power Systems (3),2415–2428 May (2018).[12] Gomez-Herrera, J. A. and Anjos, M. F.Optimization-based estimation of power capacityprofiles for activity-based residential loads. International Journal of Electrical Power & EnergySystems , 664–672 jan (2019).7UROGEN 2019 September 12-14, 2019, Guimarães, Portugal N u m b e r o f s o l u t i o n s Number of feasible solutions with c > 0 at different times Figure 10: Number of non-zero capacity solutions[13] Koliou, E., Eid, C., Chaves-Ávila, J. P., and Hakvoort,R. A. Demand response in liberalized electricitymarkets: Analysis of aggregated load participation inthe german balancing mechanism.
Energy , 245 –254 (2014).[14] Dheeru, D. and Karra Taniskidou, E. UCI machinelearning repository. http://archive.ics.uci.edu/ml , (2017).[15] Hébrail, G. and Bérard, A. Individualhousehold electric power consumption data set. https://archive.ics.uci.edu/ml/datasets/individual+household+electric+power+consumption , 08 (2012).[16] Bezanson, J., Edelman, A., Karpinski, S., and Shah,V. Julia: A fresh approach to numerical computing. SIAM Review (1), 65–98 (2017).[17] Hunter, J. D. Matplotlib: A 2d graphics environment. Computing In Science & Engineering (3), 90–95(2007).[18] JuliaStats. KernelDensity.jl v0.4.1. https://github.com/JuliaStats/KernelDensity.jl , 02(2018).[19] Lougee-Heimer, R. The Common OptimizationINterface for Operations Research: PromotingOpen-source Software in the Operations ResearchCommunity. IBM J. Res. Dev. (1), 57–66 January(2003).[20] Dunning, I., Huchette, J., and Lubin, M. JuMP: AModeling Language for Mathematical Optimization. SIAM Review (2), 295–320 (2017). [21] Riseth, A. N., Besançon, M., and Naiem, A.MultiJuMP.jl. https://doi.org/10.5281/zenodo.1042586 , 02 (2019).[22] JuliaCI. BenchmarkTools.jl v0.3. https://github.com/JuliaCI/BenchmarkTools.jlhttps://github.com/JuliaCI/BenchmarkTools.jl