Successive Convexification Algorithms for Optimizing Power Systems With Energy Storage Models

Energy storage systems (ESSs) are increasingly used in power system optimization. Different ESS mathematical models are developed that consider nonlinear functions for power losses. However, these models require non-convex constraints to represent the ESS losses, resulting in challenging optimization problems. To reduce the complexity, convex relaxation models are often derived but generate infeasible solutions when the relaxation exactness is violated. To deal with this issue, this work develops two successive convexification algorithms that generate fast and high-quality feasible solutions when the derived solution is not exact. The first algorithm handles general loss functions, while the second algorithm enhances performance when piecewise-linear loss functions are used. Specifically, the algorithms reduce the feasible region of the relaxed ESS models using a tightening box trust region around the current solution in successive iterations. The proposed algorithms are applied to the Unit Commitment and Peak Shaving and Energy Arbitrage problems to investigate their performance considering piecewise-linear and quadratic ESS loss functions. Simulation results demonstrate the impact of the ESSs relaxation violation on the actual system operation and validate the algorithms efficacy to generate high-quality feasible and even optimal solutions with significantly lower execution times compared to problems utilizing exact ESS models.

Abstract-Energy storage systems (ESSs) are increasingly used in power system optimization.Different ESS mathematical models are developed that consider nonlinear functions for power losses.However, these models require non-convex constraints to represent the ESS losses, resulting in challenging optimization problems.To reduce the complexity, convex relaxation models are often derived but generate infeasible solutions when the relaxation exactness is violated.To deal with this issue, this work develops two successive convexification algorithms that generate fast and high-quality feasible solutions when the derived solution is not exact.The first algorithm handles general loss functions, while the second algorithm enhances performance when piecewise-linear loss functions are used.Specifically, the algorithms reduce the feasible region of the relaxed ESS models using a tightening box trust region around the current solution in successive iterations.The proposed algorithms are applied to the Unit Commitment and Peak Shaving and Energy Arbitrage problems to investigate their performance considering piecewise-linear and quadratic ESS loss functions.Simulation results demonstrate the impact of the ESSs relaxation violation on the actual system operation and validate the algorithms efficacy to generate high-quality feasible and even optimal solutions with significantly lower execution times compared to problems utilizing exact ESS models.
Index Terms-Convex relaxation, energy storage models, optimization, peak shaving and energy arbitrage, unit commitment.

I. INTRODUCTION
E NERGY storage systems (ESSs) is an emerging tech- nology able to compensate the negative effects of the high penetration of renewable energy sources (RES) into the The authors are with the KIOS Research and Innovation Center of Excellence and the Department of Electrical and Computer Engineering, University of Cyprus, 1678 Nicosia, Cyprus (e-mail: ltziov01@ucy.ac.cy; lhadji02@ucy.ac.cy; stimo@ucy.ac.cy).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TSG.2023.3316720.
Digital Object Identifier 10.1109/TSG.2023.3316720power system.In general, ESSs can provide several functionalities in the electricity sector including the provision of grid services, e.g., frequency control and peak shaving, as well as the management of energy resources for electricity cost reduction [1], [2].Since these ESS functionalities allow for increased and effective RES integration, the total electricity storage capacity is expected to rapidly grow to meet the target of the European Union to become climate-neutral by 2050 [3].
The increasing utilization of ESSs in power system applications necessitates the use of mathematical models for their representation in power system optimization problems.Different ESS models have been proposed that utilize non-convex constraints to represent the ESS power losses, e.g., using piecewise linear [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18] and quadratic [4] loss functions, resulting in non-convex optimization problems.To solve these challenging problems, relaxed ESS models are often derived by relaxing the non-convex constraints.The relaxed ESS models reduce the problem complexity because they are convex, making them suitable for real-time online energy management, e.g., using model predictive control.However, the relaxed models yield infeasible solutions when the derived solutions are not exact.
The relaxed ESS models are also used in non-convex mixedinteger optimization problems, which are hard to solve, to reduce their computational complexity by decreasing the number of binary variables [12], [13], [14], [15], [16].Specifically, the relaxed ESS models are utilized in various unit commitment (UC) problems that consider electric vehicles [12], stochasticity of RES [13], and contingencies [14], as well as in expansion planning (EP) problems [15], [16].Although the relaxation exactness is shown to hold under some conditions in [14] and [15], two examples where the relaxation exactness is violated in the UC and EP problems are presented in [17].Two approaches are presented in the literature for dealing with potential violation of the exactness of the ESS relaxed constraints.The first approach is to prove their exactness under all or specific conditions to guarantee feasibility of the considered problems.The relaxation exactness is ensured for different formulations under some sufficient conditions [8], [11], [12], [18], but these conditions limit the applicability of the relaxed ESS models.The second approach is to develop methodologies that find feasible solutions when the derived solutions result in non-exact ESS constraints.In [9], a penalty term is added in the optimization objective function to avoid simultaneous charging and discharging, generating feasible solutions.However, this approach is applicable only to a specific piecewise-linear ESS model and requires a proper tuning of the penalty parameter because small values may yield infeasible solutions, while large values may lead to non-optimal solutions.In fact, tuning needs to be performed separately for each problem instance, limiting the applicability of this method.This work aims to develop successive convexification algorithms to find high-quality feasible solutions considering general ESS loss functions.
Successive approximation methodologies are widely used to solve power system optimization problems with nonlinear and non-convex constraints.Specifically, successive linear programming (SLP) approaches are utilized to solve the optimal power flow [19], [20], [21], [22], [23] and optimal power and gas flow [24] problems.The SLP is an iterative approach that linearizes the nonlinear constraints around the current feasible solution, i.e., using first-order Taylor series expansions, and solves the approximate linear program in each iteration.Sequential quadratic programming (SQP) is a successive approximation methodology similar to SLP and is also used to solve the power flow problem [25], [26], [27].Although the SLP and SQP methods are effective for solving large nonlinear optimization problems, the performance of these methods in terms of solution quality and execution time depends on the initialization of the SLP [23] and SQP [25] algorithms.To improve the computational performance of the SLP algorithm, various techniques are presented in [28], [29] to initialize the algorithm with a near-optimal solution.
This work develops two convexification algorithms that yield fast and high-quality feasible solutions when the derived solution using relaxed ESS models is not exact.The first algorithm considers various ESS models with different loss functions, while the second specialized algorithm enhances the algorithm performance when piecewise-linear ESS models are used.Specifically, both algorithms reduce the feasible region of the relaxed ESS models in each iteration using a trust region around the obtained solution.The performance of the proposed algorithms is investigated by considering two optimization problems in power systems that incorporate ESSs: (a) the Unit Commitment and (b) the Peak Shaving and Energy Arbitrage (PSEA).These problems were appropriately selected to study the impact of the relaxation violation on the actual system operation, because they can generate non-exact solutions when the relaxed ESS models are used.Simulation results validate the effectiveness of the proposed algorithms to yield high-quality feasible and even optimal solutions with significantly lower execution times compared to problems utilizing exact ESS models.In addition, the performance of the proposed algorithms is compared with literature-based approaches, i.e., the method proposed in [9] and the SQP.To the authors best knowledge, this is the first work that develops convexification algorithms that deal with potential violation of the exactness of the ESS relaxed constraints considering general power loss functions.In summary, this work has three main contributions.
1) Development of a general successive convexification algorithm that yields fast and high-quality feasible solutions considering general power loss functions.2) Development of a specialized successive convexification algorithm that enhances the solution quality and execution speed for piecewise-linear loss functions.Specifically, this algorithm is an extension of the general successive convexification algorithm to enhance the algorithm performance when piecewise-linear ESS models are used.3) Application of the proposed algorithms in two different optimization problems in power systems to investigate their performance considering piecewise linear and quadratic ESS models.The rest of the paper is organized as follows.Section II states the examined problem and Section III describes the considered ESSs models.The proposed algorithms are presented in Section IV and the UC and PSEA problems are formulated in Section V. Simulation results are shown in Section VI and conclusions are given in Section VII.

A. Generic ESS Model
The generic ESS model used in power system optimization problems is expressed in discrete time and describes the evolution of the ESS state-of-charge (SoC) 1 over time based on the charging/discharging power and power losses, given by where T = {1, . . ., T} denotes the considered time horizon and T the time-slot length in hours; K = {1, . . ., K} denotes the set of ESSs, and K the number of the considered ESSs.Variables C t,k , P S t,k and P L t,k denote the ESS SoC, discharging (P S t,k ≥ 0) or charging (P S t,k < 0) power, and power losses, respectively.The ESS energy and power limits are set as where constants C k and C k denote the minimum and maximum SoC limits, P S t,k and P S t,k the discharging and charging power limits, and I k the initial SoC.

B. ESS Power Losses
Several ESS models have been proposed in the literature that consider different functions to represent the ESS power losses in (1), depending on the ESS technology.
Lossless ESS Model: This is the simplest but most unrealistic model given by Piecewise-linear ESS Model: This is the most widely used model that represents electrochemical storage technologies [4].It is defined as where e k is a positive power loss coefficient.The absolute value in ( 4) is used to avoid negative losses when P S t,k < 0; thus, the power losses are represented by two piecewise-linear segments.Note that the model in (4) can be reformulated using a constant charging/discharging efficiency and two separate variables for charging and discharging [5], [7], [8], [11]; however, this mathematical reformulation requires the inclusion of non-convex complementarity constraints to ensure non-simultaneous charging and discharging.A variation of the piecewise-linear model in (4) to represent a flywheel energy storage is derived in [30], where the power losses are proportional to |P S t,k | and C t,k where e c k is the losses coefficient associated with the SoC.Quadratic ESS Model: The quadratic model is also used for electrochemical storage [4], given by where e q k is a positive losses coefficient.Using (6), high charging/discharging power rates are "penalized", because increased power losses are generated.
P 2 /C ESS Model: The P 2 /C model is presented in [4], [31] and is used for storage technologies for which the open circuit voltage varies with the energy stored in the ESS [4], given by General ESS Model: To represent various power loss functions using one ESS model, we define the general model as where

C. Optimization Formulation
The general ESS model ( 1)-(2b), (8) can be used in an optimization problem to allow the optimal energy scheduling of ESSs incorporated in a power system application.Let P S , P L , and C denote the vector forms of the variables P S t,k , P L t,k , and C t,k , ∀t ∈ T , k ∈ K, respectively.Let also vector y denote all the ESS variables, y = {P S , P L , C}, and vector x denote the rest design variables of a problem under consideration.Such an optimization problem, can be stated as 2b), (8) (9d) where f (x, y) is the objective function and g j (x, y) and l h (x, y) are the inequality and equality constraints associated with the considered problem, respectively.Problem ( 9) is convex, which can be fast and reliably solved, when (a) f (x, y) and g j (x, y) are convex functions and (b) l h (x, y) and ĝ(P S t,k , C t,k ) are affine functions [32].However, constraint ( 8) is non-convex in the general case because function ĝ(P S t,k , C t,k ) is nonlinear; hence, Problem ( 9) is non-convex and hard to solve.To address this issue, (8) can be relaxed to an inequality constraint, yielding the relaxed problem ) is convex, e.g., using the power loss functions of the piecewise linear, quadratic, and P 2 /C ESS models.Therefore, the relaxed problem deals with power losses using nonlinear convex constraints.The relaxed problem generates the optimal solution of the exact problem (9) when the relaxation exactness is satisfied, i.e., equality is attained in constraint (10c).Otherwise, the solution is infeasible because increased ESS power losses are generated, which means that more energy is wasted than prescribed by the power loss function.In addition to generating fast and reliable solutions when the relaxed problem is convex, a significant computational time reduction can also be achieved even when the formulated relaxed problem is a mixed-integer program (MIP), i.e., when x includes integer variables.
Obtaining equality in constraint (10c) when solving the relaxed problem (10) implies relaxation exactness and hence optimality for problem (9).Nonetheless, little attention has been given to cases where the relaxation is not exact.This work aims to develop solution methodologies that yield fast and high-quality solutions to (9) under different ESS models, even when the solution of ( 10) is infeasible for (9).
Remark: Problem (9) can include a relaxed AC optimal power flow model, e.g., the second-order cone programming (SOCP) power flow model [5], to consider network constraints in power system applications.However, a different methodology should be used to deal with potential violation of the relaxed power flow constraints.

III. ENERGY STORAGE MODELS
This section formulates the exact and relaxed versions of the piecewise-linear and quadratic ESS models and presents the feasible region of their power losses.

A. Piecewise-Linear ESS Models
Exact Model E L : The loss function in ( 4) is reformulated to define the discharging and charging power losses [5 where e d k and e c k denote the positive discharging and charging losses coefficients.Note that P L,d t,k ≥ 0 and P L,c t,k ≤ 0 when P S t,k ≥ 0 (discharging), while P L,d t,k < 0 and P L,c t,k ≥ 0 when P S t,k ≤ 0 (charging).Therefore, the power losses are defined as the maximum between P L,d t,k and P L,c t,k , given by Constraint ( 12) is non-convex and logical constraints with binary variables are needed to represent it.An alternative way to handle constraint ( 12) is by replacing P S t,k with separate variables for the charging and discharging power, P c t,k ≥ 0 and P d t,k ≥ 0, defined as Using ( 13), the power losses in ( 11)-( 12) are reformulated as The non-convex complementarity constraint (14b) ensures non-simultaneous charging and discharging, such that P L t,k = e d k P d t,k and P L t,k = e c k P c t,k when P S t,k ≥ 0 and P S t,k ≤ 0, respectively.Replacing ( 13)-(14a) in (1) yields the widelyused model formulation [7], [8], [11], for all t ∈ T , k ∈ K where constants k denote the discharging and charging efficiency coefficients.The complementarity constraints (14b) can be modelled using type 1 special ordered set (SOS-1)2 constraints, as SOS-1(P d t,k , P c t,k ).Relaxed Model R L : Deriving the convex hull of constraint (12) [5] yields the relaxed model as, for all t ∈ T , k ∈ K where constant Affine constraints (16b) and (16c) provide lower and upper bounds on the power losses according to (12), creating the feasible region shown in Fig. 1(a).Similarly, the red dashed lines provide upper bounds on the power losses, such that the shaded areas are the feasible regions of the relaxed models.For ease of representation we omit indices t, k from the variables in the figure.
Approximate Model E X : To reduce the computational complexity of Model E Q , the non-convex quadratic constraint (6) is approximated using a piecewise-linear function with N linear segments.Specifically, N + 1 points, N = {1, . . ., N + 1}, on the Cartesian plane with coordinates (x n , ŷn ), n ∈ N are selected, where xn and ŷn correspond to the values of P S t,k and P L t,k , respectively.Since two adjacent points can be used as the endpoints that represent a linear segment, the approximate model is given by Constraints (1)−(2b) (18a) where λ t,k,n are positive continuous variables which correspond to one specific point.The type 2 special ordered set (SOS-2)3 constraints (18e) ensure that the points lie on the piecewise linear curve [33].Although approximating the quadratic loss function selecting a large number of linear segments reduces the approximation error, the computational complexity increases.Fig. 2 shows an example of the piecewise-linear approximation with N = 4 and N = 20.
Relaxed Model R Q : Deriving the convex hull of the feasible set of constraint ( 6) yields the following relaxed model, for all t ∈ T , k ∈ K, where constant αt,k = (e The convex quadratic and affine constraints (19b) and (19c) provide lower and upper bounds on the power losses according to (6), creating the feasible region of Fig. 1(b).

C. Relaxation Exactness and Tightness
The Relaxed Model R L provides the same solution with the Exact Model E L when equality is attained in one of the two constraints in (16b), otherwise increased power losses occur.Similarly, the Relaxed Model R Q is exact when the equality is attained in (19b).Fig. 1 depicts the feasible region of the power losses for both the exact and relaxed versions of the piecewise linear and quadratic ESS models.As can be observed, Models R L and R Q are tight because the feasible region of these models is the convex hull of the feasible region of the corresponding Exact Models E L and E Q .

IV. SOLUTION METHODOLOGY
This section proposes a methodology to generate fast and high-quality feasible solutions, using the relaxed problem (10), when the relaxed ESS constraint (10c) is not exact.Specifically, two successive convexification algorithms for (a) general ESS models (SCA-GN) and (b) piecewise-linear ESS models (SCA-PL) are developed.Algorithm SCA-PL is a specialized version of SCA-GN to enhance the algorithm performance when piecewise-linear loss functions are used.Thus, only Model R L and both Models R L and R Q can be used in Algorithms SCA-PL and SCA-GN, respectively.Towards this direction, we define the optimization formulations P R L and P R Q that integrate Models R L and R Q , respectively, Algorithm SCA-GN: The main idea of Algorithm SCA-GN is to reduce the relaxed feasible region of the ESS power losses, depicted in Figs.1(a) and (b), in successive iterations using a tightening box trust region around the current solution.Specifically, as shown in Figs.3(a)-(c), the relaxed feasible region is reduced by adjusting the maximum discharging and charging power limits in iteration q, P S(q) t,k and P S(q) t,k , based on a trust region around the ESS power set-points, P S(q−1) t,k , obtained in the previous iteration q−1.Feasible solutions are obtained for Algorithm SCA-GN when the relaxation exactness condition is satisfied ∀t ∈ T , k ∈ K, defined as where parameter ρ ≥ 0 denotes the maximum approximation error.Note that feasible solutions are generated according to (8) when ρ = 0.However, similarly with the Approximate Model E X , an approximate solution can also be Algorithm 1 SCA-GN 1: Input: ρ, σ , , i. 2: Init.q = 0. 3: Solve P R i to obtain x (q) and y (q) (P 20) is satisfied then 5: Return x = x (q) , ŷ = y (q) .6: else Set q = q + 1. 10: Set 12: 14: Set P S(q) t,k = max(P S t,k , P S(q) t,k ), ∀t ∈ T , k ∈ K.

15:
Solve P R i to obtain x (q) and y (q) .16: Return x = x (q) , ŷ = y (q) .obtained when ρ > 0. For example, an approximate solution for Algorithm SCA-GN is presented in Fig. 3(c) when ρ = 0.08.Algorithm 1 summarizes the proposed procedure utilizing Problem P R i , i = {L, Q}.Initially, Problem P R i is solved (Line 3) such that the obtained solution is feasible when condition (20) is satisfied; otherwise, the iterative procedure is executed (Lines 7-16).Algorithm SCA-GN considers a "trust region" length L (q) t,k that defines the relative distance between the maximum charging and discharging bounds in iteration q, P S(q) t,k and P S(q) t,k , initialized as t,k is reduced in each iteration q based on the coefficient 0 < σ ≤ 1 that controls the "trust region" length (Line 10).Next, the new bounds P S(q) t,k and P S(q) t,k are calculated according to Lines 11-12 based on the reduced region length L (q) t,k and the ESS power set-points, P S(q−1) t,k , obtained in the previous iteration q − 1.Specifically, the reduction of L (q) t,k in each iteration reduces the relative distance between the charging and discharging bounds, causing the progressive reduction of the relaxed feasible region of the ESS power losses.Lines 13-14 ensure that P S(q) t,k and P S(q) t,k are within the ESS power limits P S t,k and P S t,k .Considering the example presented in Figs.3(a [11][12][13][14].Finally, Problem P R i is solved using the new bounds P S(q) t,k and P S(q) t,k (Line 15).Algorithm 1 converges when condition ( 20) is satisfied or when all L (q) t,k , ∀t ∈ T , ∀k ∈ K, become smaller than the algorithm tolerance (Line 8).Note that Algorithm SCA-GN can be used for any nonlinear convex relaxed ESS model including ( 4)- (8).
The solution quality and execution speed of Algorithm SCA-GN is dependent on the coefficients ρ, and σ .Small values of ρ are desirable because they reduce the approximation error, but increase the number of iterations and hence the execution time of Algorithm SCA-GN.Note that the tolerance is only used to ensure convergence of Algorithm 1 when condition (20) cannot be satisfied, yielding approximate solutions; thus, should be selected sufficiently small.Algorithm 1 converges in a finite number of iterations because the relaxed feasible region is progressively reduced based on the coefficient σ that controls the "trust region" length.The careful selection of the value of σ is critical because it provides a trade-off between execution speed and solution quality.For example, small values of σ increase the number of iterations and thus the execution time but enhance the solution quality.Although this method does not guarantee optimality, values in the range σ ∈ [0.5, 0.7] result in a small number of iterations (2-10) and yield high quality solutions in practical applications.
Algorithm SCA-PL: Algorithm 2, which is an extension of Algorithm 1, only utilizes Problem P R L and varies the bounds P S(q) t,k and P S(q) t,k until one of the two parameters become negative or positive, respectively.Then, the corresponding parameter is set equal to zero and the other parameter is set equal to its initial value (Lines 2-5) to consider the whole charging or discharging segment, improving the algorithm performance in terms of solution quality and execution speed.An example of Algorithm SCA-PL using the Model R L is shown in Figs.3(d)-(f).
Remark: Algorithm SCA-GN can be used for any convex power loss function, ĝ(P S t,k , C t,k ), and Algorithm SCA-PL for any piecewise-linear loss function, e.g., see (5).This can be achieved by (a) relaxing the non-convex ESS power losses constraint (8) to the convex constraint (10c), and (b) introducing upper bounds on the power losses in the formulation of the relaxed problem such that Function g(P S t,k , C t,k ) must be constructed in a way such that constraints (10c) and ( 21) define the convex hull of the exact feasible region, defined by (8).

V. POWER SYSTEM OPTIMIZATION PROBLEMS
The performance of the presented ESS models and the proposed algorithms is investigated by considering two different optimization problems in power systems, the UC and PSEA.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

A. Unit Commitment
The UC problem schedules the generating resources to satisfy the load demand over a time horizon by minimizing the total cost of operation.Based on the simplified UC problem in [17], this problem schedules the conventional generating units and ESSs to ensure the power balance between generation and demand, including ramping constraints of the units.Let G = {1, . . ., G} denotes the set with the generating units and variables P G t,g and z t,g ∈ {0, 1} denote the generating power and on/off status of the unit g ∈ G at time t ∈ T , respectively.The objective is to minimize the quadratic cost functions of units.The considered problem, denoted by U is given by minimize T t∈T g∈G αg z t,g + βg P G t,g + γg P G t,g subject to: where constants αg , βg , and γg denote the coefficients of the quadratic cost function of unit g ∈ G.In objective (22a), the fixed cost αg is included in the objective only when the unit is on, i.e., z t,g = 1.Constraint (22b) ensures that the power generation of unit g ∈ G at time t ∈ T is between its minimum and maximum limits (P G g , P G g ) when z t,g = 1; otherwise, P G t,g = 0 implying that the unit is off.Ramp-up constraints in (22c) limit the power increment of unit g ∈ G between two consecutive time intervals, where constants R U g and R SU g denote the generation upward and start-up ramp rates.Similarly, rampdown constraints are set in (22d), where constants R D g and R SD g denote the generation downward and shutdown ramp rates.Constraint (22e) ensures the power balance between produced power from the units, ESSs discharging/charging power, and load demand Dt .Similar to [17], start-up and shutdown costs, minimum up and down times, and reserves are neglected.The following five optimization problems are derived by integrating the ESS models in Problem U: • U E L and U E Q : Use the Exact Models E L and E Q .
in (22e).Utilizing the Relaxed Problems U R L and U R Q in Algorithms SCA-GN and SCA-PL yield three additional problems: and U X are mixed-integer quadratic programs (MIQPs) with SOS-1, bilinear, and SOS-2 constraints, respectively.Problems U R L , U A1 L , and U A2 L with the relaxed ESS models are MIQPs, while U R Q and U A1 Q are mixed-integer quadratically constrained quadratic programs (MIQCQPs).Problems U R L and U R Q have lower computational complexity compared to Problems U E L and U E Q and are expected to yield lower execution times.Although the UC problem is an MIP which is challenging to solve, the relaxed or exact versions of nonlinear ESS models, e.g., piecewise linear ESS model, are widely incorporated in UC formulations to represent the ESS power losses [13], [35], [36].

B. Peak Shaving and Energy Arbitrage
This problem considers the integration of a single ESS, K = {1}, in a medium voltage (MV) radial distribution grid for energy arbitrage purposes and the provision of peak shaving services to a high-to-medium voltage (HV/MV) power transformer.As shown in Fig. 4, an ESS is utilized to absorb [provide] power from the distribution grid to eliminate the reverse [direct] power limit violations of the power transformer, enabling an increased RES penetration and load demand growth [30].Specifically, the ESS can regulate the power flows through the distribution lines to avoid the overloading of the transformer and power flow limit violation of grid lines.Moreover, the ESS is used to maximize profits through energy arbitrage in electricity markets by buying and storing energy when prices are low and selling when prices are high [37].Towards this direction, an optimization problem is formulated to maximize the arbitrage profit and minimize the square of the transformer violated power, x t , ∀t ∈ T , considering the DC power flow equations to ensure the power flow limits of the distribution grid.The effectiveness of different objectives to penalize the violated power is studied in [38] Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where E denotes the set of grid lines, A the set of all buses, M j ∈ A the set of buses that are connected to bus j through grid lines, D and V the set of loads and RES plants (e.g., photovoltaics), and Dj ∈ D, Vj ∈ V and Kj ∈ K the set of loads, RES plants, and ESSs located at bus j.Variables P F t denote the transformer power, P W t,ij the power flow through the line (i, j) ∈ E, and t,j the voltage angle at bus j.Constants ĉt denote the electricity price for buying and selling power at time t, P F , P F the reverse and direct power limits of the transformer, B ij the susceptance of line (i, j) ∈ E, P D t,d the load demand d ∈ D, P V t,v the power generation of the RES plant v ∈ V, and P W ij the power flow limit of line (i, j) ∈ E. Soft constraints (23b) aim to restrain the transformer power within its limits by penalizing variables x t in the objective (23a) with a penalty coefficient W. To protect the transformer, the value of W must be sufficiently large such that the arbitrage profit is maximized as far as it does not create power violations.Constraints (23c) define the power flow through line (i, j) ∈ E as the summation of the power flows through the connected lines and the power from the loads, RES plants, and ESSs which are connected at bus j.Constraints (23d) associate the power flows with the voltage angles and (23e) determine the power flow limits of grid lines.The transformer power is defined in (23f) and the voltage angles of the reference bus are set equal to 0 in (23g).The following five optimization problems are derived by incorporating the ESS models in Problem S: • S E L and S E Q : Use the Exact Models E L and E Q .

VI. SIMULATION RESULTS
This section evaluates the performance of the ESSs models and the proposed algorithms applied in the UC and PSEA problems.All problems in Section V are coded in MATLAB and solved using optimization solver Gurobi [40], on a personal computer with 16 GB RAM and an Intel Core-i7 2.11 GHz processor.The algorithm parameters ρ and are set to 0.001 and 0.0001, respectively.To examine the solution quality of Problems U R i , U A1 i , U A2 L , S R i , S A1 i , and S A2 L , i = {L, Q}, the optimality gap metric is considered that defines the relative distance between their solution value and the where the optimal value is obtained by solving the corresponding problems with the exact ESSs models.The following cases can be observed depending on the value of the optimality gap: 1) Optimality gap = 0: The optimal solution is generated if the ESSs relaxation is exact.2) Optimality gap > 0: A feasible non-optimal solution is derived.3) Optimality gap < 0: The ESSs relaxation is non-exact and hence an infeasible solution is generated.

A. Unit Commitment
Setup: The simulation setup is composed of 2 conventional units and 6 ESSs presented in Tables I and II as well as the 24-hour load demand shown in Table III, defined as Scenario S 1 .Note that the ESSs efficiencies η d k and η c k and losses coefficients e q k , presented in Table II, are used in the piecewise linear and quadratic ESS models, respectively.In Table IV we consider Scenarios S 1 − S 6 with an increasing number of units and ESSs by (a) duplicating the units and ESSs of S 1 and (b) multiplying the load demand pro-Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.file with the load magnitude of Table IV.The scenarios were selected in a way to yield non-exact solutions to Problems U R L and U R Q .ESSs Relaxation Violation: This case study examines the impact of the ESSs relaxation violation on the solution of the UC problem, considering only the piecewise-linear ESS model.Towards this direction, we consider Scenario S 1 but use the load demand given in Table V with a time horizon of 5 hours, instead of 24 hours, for better visualization of the results.The optimal results obtained from solving Problem U E L are demonstrated in Figs.5(a L is indicated in Fig. 5(h), for the interval [3,4] h, where the energy losses resulting from the Fig. 6.Total power imbalances occurring from the ESSs relaxation violation.The total power imbalances are defined as the scheduled ESSs discharging and charging power minus the actual ESSs power.optimization solution are higher than the actual losses. 4As a result, the commitment of the "expensive" Unit 2 is avoided, yielding an operating cost of 484.26 e which is 7.65% (optimality gap = −7.65%)lower than the optimal value of 524.36 e.

TABLE V LOAD DEMAND II
To study the impact of the relaxation violation on the real system operation, we calculate the actual ESSs power setpoints 5 assuming that the charging/discharging decisions (see Fig. 5(f)) are applied in actual ESSs with embedded primary controllers.Fig. 6 shows a total power imbalance of −0.27 MW in the interval [3,4] h due to the relaxation violation, because the primary controllers reduce the ESSs charging power by 0.27 MW to ensure the SoC upper limits.Therefore, the relaxation violation causes a mismatch between scheduled and actual operation of the real system, creating power imbalances in the UC problem.The power imbalance leads to infeasibility because Unit 1 cannot reduce its generation under 23.54 MW and in the next interval to produce 38.54 MW due to its ramp-up limit of 15 MW/h.Although power imbalances in real operation can be compensated by maintaining reserves, they are undesirable because they threat the safe operation of the system.
Literature-based Method: To find a feasible solution when the ESSs relaxation exactness is violated, we employ the literature-based method proposed in [9] and presented in the Appendix.Considering the relaxed version of the piecewiselinear ESS model, a penalty term μ is utilized to penalize the ESSs charging power in objective (25a), aiming to ensure nonsimultaneous charging and discharging.Figs.7(a) and 7(b) illustrate the power imbalances, calculated based on the actual system operation, and optimality gap for various values of μ for the Scenario S 1 .Note that the optimality gap is calculated based on the objective (22a), ignoring the penalty term.Figs.7(a) and 7(b) indicate that feasible solutions, i.e., solutions that satisfy the relaxation exactness, are obtained for μ ≥ 22, where the power imbalances are eliminated.However, Fig. 7(b) demonstrates that a non-optimal solution with an optimality gap of 0.48% is generated for μ = 22, while larger values of μ increase the optimality gap to values higher than 4%.The main drawback of this method is in finding the best  L ) when σ ≤ 0.9; infeasible solutions are obtained for σ = 0.999.Algorithm SCA-PL (U A2 L ) always generates feasible solutions and yields optimal solutions for σ ≤ 0.9, indicating the superiority of Algorithm SCA-PL compared to SCA-GN when the piecewise-linear ESS model is used.When the quadratic ESS model is utilized, Algorithm SCA-GN (U A1 Q ) yields (a) feasible solutions for σ ≤ 0.7 and (b) optimal solutions for σ ≤ 0.5.Table VII shows the decrement of the algorithms iterations q as the value of σ increases, indicating that optimal solutions with a small number of iterations (q ≤ 5) are obtained for 0.5 ≤ σ ≤ 0.7 (see Table VI).
Algorithms Performance: The execution times and solution quality of the proposed algorithms are investigated for Scenarios S 2 − S 6 , for five random demand instances.Considering the quadratic ESS models, Fig. 8(a) depicts the low average execution times of Problem U R Q compared to the Q exceed the maximum time-limit set (60 minutes) for all scenarios.In Problem U X , we selected N = 20 to yield solutions with a maximum approximation error of ρ = 0.001 according to (20).As shown in Fig. 8(a), Algorithm SCA-GN for σ = 0.5 and σ = 0.7 yields significantly lower execution times compared to U X , achieving a total time reduction of 56.9% and 70.9%.This significant time reduction is achieved despite the fact that the relaxed problem remains an MIP which is challenging to solve.The figure also indicates that increasing σ results in faster solutions because the algorithm requires fewer iterations.Table VIII shows the average objective values of the considered problems for S 2 − S 4 , indicating that Problem U R Q generates solutions with lower objective values compared to U X and U A1 Q due to the relaxation violation.Interestingly, Table VIII depicts that Algorithm SCA-GN for σ = 0.5 and σ = 0.7 always yields better solutions compared to U X , despite the fact that U X is almost optimal.This clearly indicate the high quality of the obtained solutions.Considering the piecewise linear ESS models, Algorithm SCA-PL for σ = 0.7 yields considerably lower execution times compared to Problem U E L , as shown in Fig. 8(b) for S 4 − S 6 , achieving a total time reduction of 49.85%.Table IX shows the average objective values of the considered problems for S 4 − S 6 , indicating that Problem U R L yields non-exact solutions with lower objective values compared to the optimal values obtained using U E L (optimality gap = -0.24%).The table also indicates that Algorithm SCA-PL for σ = 0.7 yields almost optimal solutions with a maximum optimality gap of 0.00002%.The results in both Figs.8(a) and (b) show that the proposed algorithms execute significantly faster compared to the solution of exact problems with stateof-the-art solvers when the total number of ESSs incorporated in the UC problem increases.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

B. Peak Shaving and Energy Arbitrage
Setup: The simulation setup is composed of an ESS installed at Bus 6 of the 11-bus distribution grid shown in Fig. 4. The ESS has a usable capacity of 1.2 MWh and charging/discharging power of 2.4 MW, where the efficiencies η c = η d = 0.96 (e d = 0.0416, e c = 0.04) and the losses coefficient e q = 0.02 are used in the piecewise linear and quadratic ESS models, respectively.Fig. 9(a) illustrates 17 net-load curves constructed from historical data of a real distribution grid, obtained from a substation in Larnaca, Cyprus.Specifically, the net-load curves represent the aggregate power from all the loads and RES plants of the distribution grid.As shown by the net-load curves, high reverse and direct power flows are presented during the noon and evening hours due to intense photovoltaic generation and high load demand, respectively.Since power violations do not occur in the real distribution transformer due to its large size, we consider a smaller transformer with power limits of P F = 1.55 MW and P F = −0.5 MW, respectively.The power flow limit of lines is set to P W ij = 2.5 MW, ∀(i, j) ∈ E, such that congestion of grid lines is avoided.Fig. 9(b) depicts the electricity price for buying and selling power.The horizon is set to one day with 3-minute time intervals, the solver time limit is set to 15 minutes, and parameter W is set to 10000.Since Problems S A2 L and S A1 Q are convex and can be solved fast, parameter ρ is set to 0.0001, instead of 0.001, to improve the approximation error according to (20).  considering an actual ESS with an integrated primary controller.As shown in Figs.10(e) and 10(f), the actual transformer power deviates from the scheduled power at time 14:00, presenting high power violations, due to the reduction of the ESS charging power.As shown in Figs.10(g) and 10(h), the primary controller limits the ESS charging power because the actual ESS SoC reached its maximum value before the scheduled SoC due to the reduced actual losses.As a result, the ESS relaxation violation causes a mismatch between scheduled and actual operation of the real system, increasing the peak power violations, from 0.34 MW to 0.75 MW, and reducing the actual arbitrage profit, from 15.54 e to 14.78 e, compared to the optimal solution.To protect the transformer, high peak power violations can be reduced by applying RES power curtailments in real operation, but this action is undesirable because it reduces the RES penetration in the power systems.
Algorithms Performance: This case study investigates the performance of the proposed algorithms for the 17 net-load curves.Considering the piecewise-linear ESS models, Table X demonstrates that Problems S R L and S A2 L yield small execution times compared to the high times of the non-convex Problem S E L .However, Problem S R L yields infeasible solutions with an average optimality gap of −37.93% due to the relaxation violation.The maximum time of S E L indicates that the considered time limit of 60-minutes is exceeded for S E L in some scenarios, yielding sub-optimal solutions.Therefore, Algorithm SCA-PL (S A2 L ) for σ = 0.5 generated slightly better solutions compared to S E L .Table X also indicates that Algorithm SCA-PL generates solutions with an average execution time of 0.039 minutes, achieving a speedup of 368 (368x) times compared to the average time (14.35  ) and Problem S X for N = 20 yield similar average objective values, presenting a difference of 0.06%.Moreover, Algorithm SCA-GN yields better solutions compared to S E Q (SQP), reducing the average objective value by 0.51%.Algorithm SCA-GN generates solutions with an average execution time of 0.051 minutes, achieving a speedup of more than 1176 (1176x) and 646 (646x) times compared to the average time of Problems S X and S E Q (SQP).

VII. CONCLUSION
This work developed two convexification algorithms that yield fast and high-quality feasible solutions when the derived solution using relaxed ESS models is not exact.The first algorithm handles general power loss functions, while the second specialized algorithm enhances performance when piecewiselinear loss functions are used.The proposed algorithms are applied to the UC and PSEA problems considering piecewiselinear and quadratic loss functions.Simulation results show that the relaxed ESS models generate non-exact solutions that yield inaccurate operating costs and cause a mismatch between scheduled and actual operation of the real system, creating power imbalances and increasing the peak power violations in the UC and PSEA problems, respectively.Simulation results also indicate the capability of the proposed algorithms to yield almost optimal, if not optimal, solutions with significantly lower execution times compared to state-of-the-art solvers that utilize exact ESS models.Specifically, the proposed algorithms reduce the average execution time by 50% for the UC problem, which remains nonconvex even upon relaxation, and achieve a 2-3 orders of magnitude speedup for the PSEA problem, which becomes convex upon relaxation.

APPENDIX LITERATURE-BASED METHOD
Considering the relaxed version of the piecewise-linear model, the literature-based method proposed in [9] can be used to find feasible solutions when the relaxation exactness is violated.Specifically, this method adds a penalty term in the

Manuscript received 23
February 2023; revised 19 July 2023; accepted 2 September 2023.Date of publication 18 September 2023; date of current version 21 February 2024.This work is funded by the Horizon Europe Research and Innovation Programme under grant agreement No 101075747 (TRANSIT), the Cyprus Research and Innovation Foundation under Project CULTURE/ AWARD-YR/0322B/0003 (INVERGE), and supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No 739551 (KIOS CoE) and from the Republic of Cyprus through the Deputy Ministry of Research, Innovation and Digital Policy.(Corresponding author: Stelios Timotheou.)

Fig. 1 .
Fig. 1.Power losses of the exact and relaxed versions of the (a) piecewise linear and (b) quadratic ESS models as a function of the charging/discharging power.The blue solid lines have the dual role of (a) presenting the feasible region of the exact models, indicating the actual losses, and (b) providing lower bounds on the power losses with respect to the relaxed models.Similarly, the red dashed lines provide upper bounds on the power losses, such that the shaded areas are the feasible regions of the relaxed models.For ease of representation we omit indices t, k from the variables in the figure.

Fig. 2 .
Fig. 2. ESSs power losses as a function of the charging/discharging power for Models E Q and E X : (a) N = 4 and (b) N = 20 linear segments are considered for the piecewise linear approximation in Model E X .

Fig. 5 .
Fig. 5. Optimization results obtained by solving Problems U E L (a)-(d) and U R L (e)-(h): (a) and (e) generating power of the conventional units, (b) and (f) total ESSs discharging/charging power, (c) and (g) total ESSs state-of-charge, and (d) and (h) total ESSs energy losses.
)-5(d), while the results generated from solving Problem U R L are shown in Figs.5(e)-5(h).The generating power of the two units is illustrated in Figs.5(a) and 5(e), indicating that U E L commits both units at 4 − 5 h while U R L commits only Unit 1.The total ESSs discharging/charging power is presented in Figs.5(b) and 5(f) and the total ESSs SoC is depicted in Figs.5(c) and 5(g).The total ESSs energy losses based on the discharging/charging power decisions are shown in Figs.5(d) and 5(h).The ESSs relaxation violation in U R

Fig. 8 .
Fig. 8. Average execution times of the problems with the relaxed and exact ESS models as well as the proposed algorithms when the (a) quadratic and (b) piecewise linear ESS models are considered.

Fig. 9 .
Fig. 9. Input data: (a) net-load curves of the distribution grid and (b) energy price in e/MWh.The black load curve is used to study the ESS relaxation violation, while all curves are used to examine the algorithms performance.

Fig. 10 .
Fig. 10.Optimization results obtained by solving Problems S E L and S R L (a)-(d) and actual results of Problem S R L considering the ESSs operating limits (e)-(h): (a) and (e) transformer power, (b) and (f) discharging/charging power, (c) and (g) ESS state-of-charge, and (d) and (h) ESS energy losses.

TABLE IX AVERAGE
OBJECTIVE VALUES OF THE PROBLEMS THAT CONSIDER THE PIECEWISE LINEAR ESS MODELS minutes) of Problem S E L .Considering the quadratic ESS models, we solve the nonconvex Problem S E Q with the exact ESS model by employing the nonlinear optimization solver Artelys Knitro [40] and using the SQP algorithm.Towards this direction, we define Problem S E Q (SQP) as • S E Q (SQP): Solves Problem S E Q using the SQP algorithm.Table XI depicts the small execution times of the convex Problems S R Q and S A1 Q compared to the high times of the non-convex Problems S X and S E Q (SQP); however, Problem S R Q generates non-exact solutions.As shown in Table XI, Algorithm SCA-GN (S A1 Q Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.optimizationobjectivefunction to penalize the ESSs charging power, aiming to ensure non-simultaneous charging and discharging.Applying this method in the UC problem, the optimization problem of Section V is reformulated as denotes the penalty term of the charging power.The Relaxed Model R L is considered in (22f) and constraint (13) is used in Problem 25 to represent the ESS charging and discharging power using two separate variables.Note that small values of μ may yield non-exact solutions, while large values may generate non-optimal solutions.