Towards multiobjective optimization and control of smart grids

TL;DR

This paper explores multiobjective optimization for smart grids, balancing load shaping and flexibility using Pareto optimality to analyze trade-offs and improve grid robustness and resilience.

Contribution

It introduces a multiobjective framework employing Pareto optimality to analyze and manage trade-offs in smart grid energy storage control.

Findings

Pareto frontier analysis quantifies trade-offs between objectives.

Multiobjective approach improves balancing of load shaping and flexibility.

Framework enhances robustness and resilience of smart grids.

Abstract

The rapid uptake of renewable energy sources in the electricity grid leads to a demand in load shaping and flexibility. Energy storage devices such as batteries are a key element to provide solutions to these tasks. However, typically a trade-off between the performance related goal of load shaping and the objective of having flexibility in store for auxiliary services, which is for example linked to robustness and resilience of the grid, can be observed. We propose to make use of the concept of Pareto optimality in order to resolve this issue in a multiobjective framework. In particular, we analyse the Pareto frontier and quantify the trade-off between the non-aligned objectives to properly balance them.

Equations

$$\displaystyle\begin{split}x_{i}(k+1)&=&\alpha_{i}x_{i}(k)+T(\beta_{i}u_{i}^{+}(k)+u_{i}^{-}(k))\\ z_{i}(k)&=&w_{i}(k)+u_{i}^{+}(k)+\gamma_{i}u_{i}^{-}(k)\end{split}$$

$$0\le$$

$${\underline{u}}_i\le$$

$$0\le$$

$$x(k+1)=\left(\begin{array}{c}{x}_1(k+1)\\ x_2(k+1)\\ ⋮\\ x_{\mathcal{I}}(k+1)\end{array}\right)=\left(\begin{array}{c}{f}_1(x_1(k),u_1(k))\\ f_2(x{(}_2(k),u_2(k))\\ ⋮\\ f_{\mathcal{I}}(x_{\mathcal{I}}(k),u_{\mathcal{I}}(k))\end{array}\right)=f(x(k),u(k)),$$

$$w(n)=(w_1(n),w_2(n),\dots ,w_{\mathcal{I}}(n){)}^{\top },n=0,1,2,\dots ,k,$$

$$z(k)=z(k;(u(n){)}_{n=0}^k,(w(n){)}_{n=0}^k,x^0).$$

$$\bar{\zeta }(n)=\frac{1}{\mathcal{I}\cdot \min \{N,n+1\}}\sum _{j=n-\min \{n,N-1\}}^n\sum _{i=1}^{\mathcal{I}}{w}_i(j).$$

$$\mathbf{u}=(u(n){)}_{n=k}^{k+N-1}\mapsto \frac{1}{N}\sum _{n=k}^{k+N-1}{\left(\frac{1}{\mathcal{I}}\sum _{i=1}^{\mathcal{I}}\left[w_i(n)+u_i^{+}(n)+{\gamma }_i{u}_i^{-}(n)\right]-\bar{\zeta }(n)\right)}^2,$$

$$\displaystyle\underaccent{\bar}{c}(n)\quad\leq\quad\frac{1}{\mathcal{I}}\sum_{i=1}^{\mathcal{I}}z_{i}(n)\quad\leq\quad\bar{c}(n)\qquad\forall\,n\in[k:k+N-1]$$

$$\displaystyle\mathbf{u}\quad\mapsto\quad\sum_{n=k}^{k+N-1}\left[\left(\max\Big{\{}0,\frac{1}{\mathcal{I}}\sum_{i=1}^{\mathcal{I}}z_{i}(n)-\bar{c}(n)\Big{\}}\right)^{2}+\left(\max\Big{\{}0,\underaccent{\bar}{c}(n)-\frac{1}{\mathcal{I}}\sum_{i=1}^{\mathcal{I}}z_{i}(n)\Big{\}}\right)^{2}\right],$$

$$\underset{\mathbf{u}\in {\mathbb{U}}^N}{\min }$$

$$\underset{x\in \mathcal{M}}{\min }f(x)$$

$$f_i(x)<f_i(x^{\star })⟹\exists j\in \{1,\dots ,m\}:f_j(x^{\star })<f_j(x),$$

$$\underline{\mathbf{s}}=(\underline{s}(k),\underline{s}(k+1),\dots ,\underline{s}(k+N-1){)}^{\top }\text{ and }\bar{\mathbf{s}}=(\bar{s}(k),\bar{s}(k+1),\dots ,\bar{s}(k+N-1){)}^{\top },$$

$$\displaystyle\underaccent{\bar}{c}(n)-\underaccent{\bar}{s}(n)\quad\leq\quad\frac{1}{\mathcal{I}}\sum_{i=1}^{\mathcal{I}}z_{i}(n)\quad\leq\quad\bar{c}(n)+\bar{s}(n)\qquad\text{ for all $n\in[k:k+N-1]$}$$

$$\mathbf{s}\mapsto {\Vert \mathbf{s}\Vert }_2^2=\sum _{n=k}^{k+N-1}\underline{s}(n{)}^2+\bar{s}(n{)}^2.$$

$$\bar{z}(n)=\frac{1}{\mathcal{I}}\sum _{i=1}^{\mathcal{I}}{z}_i(n)\text{ for all n∈[k:k+N−1]}.$$

$$\bar{\mathbf{z}}\mapsto \frac{1}{N}\sum _{n=k}^{k+N-1}(\bar{z}(n)-\bar{\zeta }(n){)}^2=\frac{1}{N}{\Vert \bar{\mathbf{z}}-\bar{\zeta }\Vert }_2^2.$$

$${\mathbb{D}}_i$$

$$:=\left\{{\mathbf{z}}_i\in {\mathbb{R}}^N|\exists \mathbf{u}\in {\mathbb{U}}^N:\begin{array}{c}{x}_i(n+1)={\alpha }_i{x}_i(n)+T({\beta }_i{u}_i^{+}(n)+u_i^{-}(n))\in [0,C_i],\\ z_i(n)=w_i(n)+u_i^{+}(n)+{\gamma }_i{u}_i^{-}(n),n\in [k:k+N-1]\end{array}\right\}$$

$$\underset{(\bar{\mathbf{z}},\mathbf{s})}{\min }$$

$$\mathrm{subject}\mathrm{to}$$

$$\bar{\mathbf{z}}:\exists ({\mathbf{z}}_1^{\top },\dots ,{\mathbf{z}}_{\mathcal{I}}^{\top }{)}^{\top }\in \mathbb{D}\text{ satisfying eqrefNotationZbar,}$$

$$\underset{(\bar{\mathbf{z}},\mathbf{s})}{\min }$$

$$\mathrm{subject}\mathrm{to}$$

$$\bar{\mathbf{z}}:\exists ({\mathbf{z}}_1^{\top },\dots ,{\mathbf{z}}_{\mathcal{I}}^{\top }{)}^{\top }\in \mathbb{D}\text{ satisfying eqrefNotationZbar}$$

$$\underset{(\mathbf{z},\mathbf{s})}{\min }$$

$$\mathrm{subject}\mathrm{to}$$

$$\bar{\mathbf{z}}=\frac{1}{\mathcal{I}}\sum _{i=1}^{\mathcal{I}}{\mathbf{z}}_i$$

$${\mathbf{z}}_i\in {\mathbb{D}}_i\forall i\in [1:\mathcal{I}]$$

$${\mathbf{z}}_i^{\ell +1}\in \underset{{\mathbf{z}}_i\in {\mathbb{D}}_i}{\mathrm{argmin}}\frac{\rho }{2}{\Vert {\mathbf{z}}_i-{\mathbf{z}}_i^{\ell }+{\Pi }^{\ell }\Vert }_2^2$$

$$({\bar{\mathbf{a}}}^{\ell +1},{\mathbf{s}}^{\ell +1})\in \underset{(\bar{\mathbf{a}},\mathbf{s})\in \mathbb{S}}{\mathrm{argmin}}\kappa g(\bar{\mathbf{a}})+(1-\kappa )h(\mathbf{s})+\frac{\rho \mathcal{I}}{2}{\Vert {\bar{\mathbf{z}}}^{\ell +1}-\bar{\mathbf{a}}+\frac{{\bar{\lambda }}^{\ell }}{\rho }\Vert }_2^2,$$

Full text

Towards multiobjective optimization

and control of smart grids††thanks: The authors acknowledge funding by the Federal Ministry for Education and Research (BMBF; grant 05M18SIA). Karl Worthmann is indebted to the German Research Foundation (DFG; grant WO 2056/6-1).

Philipp Sauerteig  and Karl Worthmann Institute for Mathematics, Faculty of Mathematics and Natural Sciences, Technische Universität Ilmenau, Ilmenau, Germany (philipp.sauerteig@tu-ilmenau.de).Institute for Mathematics, Faculty of Mathematics and Natural Sciences, Technische Universität Ilmenau, Ilmenau, Germany (karl.worthmann@tu-ilmenau.de).

The rapid uptake of renewable energy sources in the electricity grid leads to a demand in load shaping and flexibility. Energy storage devices such as batteries are a key element to provide solutions to these tasks. However, typically a trade-off between the performance related goal of load shaping and the objective of having flexibility in store for auxiliary services, which is for example linked to robustness and resilience of the grid, can be observed. We propose to make use of the concept of Pareto optimality in order to resolve this issue in a multiobjective framework. In particular, we analyse the Pareto frontier and quantify the trade-off between the non-aligned objectives to properly balance them.

1 Introduction

The energy transition has triggered a rapid uptake of renewable energies in the electricity grid, which may lead to severe problems in the energy supply, see, e.g. the introductory paragraphs of [1, 2] and the references therein. In this paper, we focus on the integration of energy storage devices [3], e.g. batteries, in microgrids [4], see also [5] for a review article and [6] for questions concerning AC/DC. Based on forecasts [7, 8], consumption and production peaks can be anticipated, such that a receding horizon strategy, typically realized via a Model Predictive Control (MPC) scheme [9, 10, 11], is idealy suited to tackle this task, see also the textbooks [12, 13] for an introduction to (nonlinear) MPC. An aspect of particular interest is whether the optimization based control algorithm is realized in a decentralized, distributed, or centralized fashion, see, e.g. [2] and [14].

A major concern from a grid operator’s perspective is the volatile power demand due to residential energy generation. In recent papers, the authors exploit inherent flexibility provided by energy storage devices to smoothen the aggregated power demand profile by approximating a given reference value [2]. Here, a key issue was the realization of the global optimum by using distributed control — both in a cooperative [15] and a non-cooperative setting [16]. Another – potentially conflicting – objective is to stay within time-varying tubes as introduced in [17], which may not contain the reference value enforced by the first optimization goal. The motivation behind this approach is to provide flexibility to the grid operator, which may be used to counteract other shortcomings within the network [18]. The contribution of this paper is to embed this problem into a multiobjective framework, link it to the scalarized objective function and investigate the connection to the Pareto frontier. Furthermore, we also introduce a more restrictive optimality concept – proper optimality in the sense of Geoffrion. By doing so, it becomes possible to quantify the trade-off between the different objectives, which is an essential prerequisite for structured decision making process.

Moreover, we provide numerical simulations to investigate the MPC closed loop within our multiobjective framework. Hereby, we use data from an Australian distribution network dataset, see [19] for further details.

Notation: For two integers $m$ and $n$ with $m\leq n$, we use the notation $[m:n]:=\{m,m+1,\ldots,n\}$.

2 Problem Formulation

In this section, we firstly introduce a mathematical model of a prosumer, i.e. a residential energy system equipped with both solar-photovoltaic panel (or any other device for energy generation) and an energy storage device like a battery. Secondly, we consider a (potentially) large network of $\mathcal{I}$, $\mathcal{I}\in\mathbb{N}$, such subsystems. The individual subsystems are connected via a common point of coupling, the grid operator (simply termed Central Entity; CE), see Figure 1 (left). So far, we essentially follow the modelling approach presented in [1, 2] and further developed in [17]. An alternative would be to also model the underlying network structure as graph and, then, to incorporate the physical connections within the microgrid, see, e.g. [20].

The CE may pursue different objectives, which leads to a Multiobjective Optimization Problem (MOP) since we assume that the subsystems are cooperative. Here, we exemplarily present two criteria and numerically investigate the network performance w.r.t. each of the two objectives.

2.1 Modelling of the System Dynamics and the Constraints

The $i$-th subsystem, $i\in[1:\mathcal{I}]$, is given by

[TABLE]

where $x_{i}(k)$ describes the State Of Charge (SOC) of the battery and $z_{i}(k)$ the power demand at time instant $k\in\mathbb{N}_{0}$. The parameter $T>0$ is the length of one time step in hours, e.g. $T=0.5$ corresponds to a half hour window. The variables $u_{i}^{+}(k)$ and $u_{i}^{-}(k)$ represent the charging and discharging rate, resp., while $w_{i}(k)$ is the net consumption (load minus generation). The constants $\alpha_{i},\beta_{i},\gamma_{i}\in(0,1]$ model efficiencies w.r.t. to self discharge and energy conversion, resp. The initial SOC at $k=0$ is denoted by $x_{i}^{0}$, i.e. $x_{i}(0)=x_{i}^{0}$ holds. The system variables are subject to the constraints

[TABLE]

where $C_{i}\geq 0$ denotes the battery capacity . Constraint (2d) ensures that the bounds on charging (2c) and discharging (2b) hold, even if charging and discharging occur during the same time interval, see Figure 1 (right). For a concise notation, we define the state constraint $\mathbb{X}_{i}:=[0,C_{i}]$ and the set of feasible control values $\mathbb{U}_{i}:=\Set{(u_{i}^{-},u_{i}^{+})^{\top}\in\mathbb{R}^{2}}{\text{{\eqref{eq:constraints:discharge}-\eqref{eq:constraints:chargeANDdischarge} hold}}}$. Furthermore, we use the notation $\bar{\cdot}$ to denote the average value over all subsystems, e.g. $\bar{w}(k)=\frac{1}{\mathcal{I}}\sum_{i=1}^{\mathcal{I}}w_{i}(k)$.

2.2 First Objective: Peak Shaving

The Central Entity is located at the common point of coupling, where the aggregated sum $\sum_{i=1}^{\mathcal{I}}z_{i}(k)$ is collected. A positive sign corresponds to a power demand while a negative sign indicates a surplus, which reflects that the aggregated local generation exceeds the aggregated local load at time instant $k$. Hence, the CE requires balancing energy to match the aggregated power demand. While that is comparatively easy for a constant power demand/supply, it is significantly more demanding if the power demand exhibits large fluctuations. Hence, a typical goal is to flatten the energy demand of the network in consideration.

For a mathematical formulation, we first compose an overall system by setting up the system dynamics

[TABLE]

the constrains $x(k)\in\mathbb{X}:=\mathbb{X}_{1}\times\ldots\times\mathbb{X}_{\mathcal{I}}\subset\mathbb{R}^{\mathcal{I}}$ and $u(k)\in\mathbb{U}:=\mathbb{U}_{1}\times\ldots\times\mathbb{U}_{\mathcal{I}}\subset\mathbb{R}^{2\mathcal{I}}$. Hence, for given data

[TABLE]

the power demand at time $k\geq 0$ resulting from a discharging/charging behaviour $(u(n))_{n=0}^{k}=((u^{-}(n),u^{+}(n)))_{n=0}^{k}\in\mathbb{U}^{k+1}$ and initial state of charge $x^{0}\in\mathbb{X}$ is

[TABLE]

Here, we emphasize that the output depends on the time varying load and generation data given by (3), which has to be predicted for load shaping, see, e.g. [21, 22] for load forecasting, [23, 24] and [25, 26] for predictions on the energy generation due to wind and solar power, resp. Hence, at time instant $k$, we optimize over a finite time horizon of $N$, $N\in\mathbb{N}_{\geq 2}$, time steps. Moreover, we require suitable (time varying) reference values, which are constructed based on the overall average net consumption, i.e.

[TABLE]

Note that $\bar{\zeta}(n)$ is the average over the previous $N$ time steps if the data history is sufficiently long, which is reflected by taking the minimum. Then, for given data $w(n)$, $n=\max\{0,k-(N-1)\},\ldots,k,\ldots,k+N-1$, and state $x(k)$, $x(k)\in\mathbb{R}^{\mathcal{I}}$, the CE wants to minimize the objective function $J_{1}:\mathbb{U}^{N}\to\mathbb{R}_{\geq 0}$,

[TABLE]

subject to the system dynamics (1) and the state constraint $x_{i}(n)\in\mathbb{X}_{i}$ for all $n\in[k:k+N]$ and $i\in[1:\mathcal{I}]$, which penalizes the deviation of the average power demand from the sequence of reference values given by the vector $\bar{\zeta}=(\bar{\zeta}(k),\ldots,\bar{\zeta}(k+N-1))^{\top}$. Here, we call $N$ the optimization horizon, on which it is assumed that the load and generation data can be reliably predicted.

The optimal control value of this optimization problem is unique and is attained as shown in [15]. Hence, we can define a state feedback law by using the first element $u^{\star}(k)$ of an optimal control sequence $\mathbf{u}^{\star}=(u^{\star}(n))_{n=k}^{k+N-1}$, which is then implemented at the plant. Then, based on the new state $x(k+1)$ and the new input data (in particular $w(k+N)$ and $\bar{\zeta}(k+N)$), the procedure can be repeated ad infinitum to generate a closed loop in a receding horizon fashion as explicated in Algorithm 1 solely based on $J_{1}$ (and, thus, ignoring the – potentially conflicting – second objective $J_{2}$, see, e.g. [2] or [12] for a more detailed introduction to (distributed) MPC).

Numerical results illustrating the outcome of such an MPC scheme are depicted in Figure 2. We point out that the reference values (4) cannot be traced due to the constraints (2).

2.3 Second Objective: Flexibility via Tube Constraints

Besides flattening the power demand profile, another goal of the grid operator may be to provide flexibility to the lower/higher layers of the overall system. The idea is to introduce time-varying tube constraints on the average power demand, i.e., we would like to achieve

[TABLE]

for the lower and upper bounds \underaccent{\bar}{\mathbf{c}}(n) and $\bar{\mathbf{c}}(n)$, $n\in[k:k+N-1]$. However, since the inclusion of the tube constraints (5) may render the problem infeasible, we minimize the violation of the tube constraints by using the objective function $J_{2}:\mathbb{U}^{N}\rightarrow\mathbb{R}_{\geq 0}$,

[TABLE]

subject to the system dynamics (1) and the state constraint $x_{i}(n)\in\mathbb{X}_{i}$ for all $n\in[k:k+N]$ and $i\in[1:\mathcal{I}]$. According to this construction a positive value of $J_{2}(\mathbf{u})$ reflects a violation of (5) which causes a loss of flexibility in the system.

Note that this cost function exhibits some complementarity: for a given time index $n$, strict positivity of one term implies that the other term is equal to zero, which reflects that a violation of Condition (5) can only occur on one side. Moreover, for a given time index $n$, the objective function is convex. Numerical solutions are displayed in Figure 3 to illustrate our modelling approach.

2.4 Multiobjective Optimization Problem and MPC Scheme

Combining the non-aligned (conflicting) objectives introduced in the preceeding subsections leads to the following multiobjective optimization problem

[TABLE]

For this kind of Multiobjective Optimization Problems (MOPs), see, e.g. [27], the following optimality concept is typically used.

Definition 1** (Pareto optimality, Pareto frontier).**

Consider the vector optimization problem

[TABLE]

with $f:X\to\mathbb{R}^{m}$, $\mathcal{M}\subseteq X\subseteq\mathbb{R}^{n}$. A vector $x^{\star}\in\mathcal{M}$ is called (Pareto) optimal or (Pareto) efficient for the optimization problem (6) if for all $x\in\mathcal{M}$ and $i\in\{1,\ldots,m\}$ the following holds:

[TABLE]

where $f_{i}:X\to\mathbb{R}$ denotes the $i$-th component of $f$, $i\in[1:m]$. Furthermore, a point $x^{\star}\in\mathcal{M}$ is called weakly (Pareto) optimal or weakly (Pareto) efficient of (6) if there is no $x\in\mathcal{M}$ with $f_{i}(x)<f_{i}(x^{\star})$ for all $i\in[1:m]$. The set $\Set{f(x)}{x\text{ is efficient}}$ is called Pareto frontier.

Remark 2**.**

Clearly, efficiency yields weak efficiency.

Before we proceed, we reformulate the MOP described above in order to facilitate its solution in Section 3. To this end, we begin with the objective function $J_{2}$. Here, we introduce the stacked auxiliary variable $\mathbf{s}=(\underline{\mathbf{s}}^{\top},\bar{\mathbf{s}}^{\top})^{\top}\in\mathbb{R}_{\geq 0}^{2N}$ with

[TABLE]

which was already used in [17]. Then, we may introduce the relaxed constraints

[TABLE]

and minimize the Euclidean norm of the vector $\mathbf{s}$. To this end, we introduce the cost function $h:\mathbb{R}^{2N}_{\geq 0}\to\mathbb{R}_{\geq 0}$,

[TABLE]

In summary, the additional signed (slack) variables, which are stacked in the vector $\mathbf{s}$, allow to represent the optimization objective $J_{2}$ by the strictly convex function $h$ at the expense of incorporating the constraints given by (7). Doing so ensures feasibility of the optimization problem, e.g. by not using the storage devices at all. If desired from an implementational point of view, the complementarity constraint $\underline{s}(n)\bar{s}(n)=0$, $n\in[k:k+N-1]$, may be included while both feasibility and the value of the objective function are maintained.

Furthermore, we introduce the auxiliary variable $\bar{\mathbf{z}}=(\bar{z}(k),\ldots,\bar{z}(k+N-1))^{\top}\in\mathbb{R}^{N}$ and the additional equality constraints

[TABLE]

Using $\bar{\mathbf{z}}$ and the constraint (8) allows us to replace $J_{1}$ by the function $g:\mathbb{R}^{N}\rightarrow\mathbb{R}_{\geq 0}$,

[TABLE]

Next, we want to rephrase the constraint set such that the output variables $\mathbf{z}_{i}=(z_{i}(k),z_{i}(k+1),\ldots,z_{i}(k+N-1))^{\top}$, $i\in[1:\mathcal{I}]$, can be used as an optimization variable instead of $\mathbf{u}$. To this end, we define the set $\mathbb{D}=\mathbb{D}_{1}\times\ldots\times\mathbb{D}_{\mathcal{I}}$ with

[TABLE]

with $x_{i}(k)\in[0,C_{i}]$ for all $i\in[1:\mathcal{I}]$. Note that, for all $i\in[1:\mathcal{I}]$, the set $\mathbb{D}_{i}$ and, thus, also the set $\mathbb{D}$ is convex and compact in view of the linearity of the constraints (1) and (2).

Then, plugging $\bar{z}$ into the constraint (7) yields the following MOP

[TABLE]

where the set $\mathbb{S}$ is convex. In (2.4) the optimization variables and, thus, the objectives are coupled via the constraint $\mathbb{S}$. If, in addition, also individual objectives of the residential energy systems have to be taken into account, the vectors $\mathbf{z}_{1},\ldots,\mathbf{z}_{\mathcal{I}}$ have to be used as optimization variables.

Finally, we explicate the Model Predictive Control (MPC) scheme to be applied.

In general, it is quite difficult to provide a meaningful analysis of the MPC closed loop resulting from a multiobjective problem formulation. However, adding a few additional assumptions allows to provide some guarantees if a stabilization task is considered, see, e.g. [28]. In our analysis, however, we restrict ourselves to the analysis of the static MOP to be solved in Step 2 of Algorithm 1 and provide numerical simulations to investigate the MPC closed loop.

3 Scalarization and Numerical Solution via ADMM

In the following, we consider the Scalarized Multiobjective Optimization Problem (SMOP) using the scalarization parameter $\kappa$, $\kappa\in[0,1]$. For details on the scalarization of MOPs, we refer to [29]. Overall, we obtain the following SMOP.

[TABLE]

In this section, we firstly present a numerical case study. Then, we briefly recall ADMM before we apply this optimization technique to solve (3). ADMM is an algorithm for distributed optimization of convex, large-scale systems, which allows to solve (3) even for a very large number $\mathcal{I}$ of subsystems. To this end, we exploit the reformulation of the original (2.4). We show that the proposed algorithm yields, for each $\kappa\in[0,1]$, not only the optimal value of the scalarized problem but that the corresponding values of the components given by $g$ and $h$ are also unique. Then, based on the results derived in this section, we investigate the connection of the parameter $\kappa$ to the Pareto frontier in the subsequent section.

3.1 Numerical Simulations

Increasing the value $\kappa\in[0,1]$ shifts the focus from penalizing the violation of the tube constraints to flatten the aggregated power demand profile. In Figure 4 the corresponding closed-loop performance at time instant $k=0$ is visualized. Here, we set $\mathcal{I}=10$, $T=0.5$ (time step of $30$ minutes), $N=48$ (one day prediction horizon), $x_{i}(0)=0.5$, $C_{i}=2$ (battery capacity), and $-\underline{u}_{i}=\bar{u}_{i}=0.5$ (maximal discharging/charging rate) and vary the tube constraints from $(\underline{\mathbf{c}}_{1},\bar{\mathbf{c}}_{1})=(0.2,0.4)$ to $(\underline{\mathbf{c}}_{2},\bar{\mathbf{c}}_{2})=(0.6,0.8)$ on a time window of three days. Note that the choice of the tube constraints seems inappropriate since it is not aligned with the uncontrolled power demand profile (dotted black line in Figure 4). However, doing so, allows us to ensure that the objectives are potentially conflicting goals. Since weather forecasting is not topic of this paper we use real world data on load and generation provided by an Australian distribution network for the prediction of the future net consumption $w_{i}$, $i\in[1:\mathcal{I}]$. The choice of $T$ corresponds to this data. In [17] the authors chose a different cost function $g$ and simply summed the objectives, which complies with a weighting parameter $\kappa=0.5$.

The corresponding open-loop perfomance, i.e. implementing $u^{\star}=(u^{\star}(k)^{\top},\ldots,u^{\star}(k+N-1)^{\top})^{\top}$ instead of $u^{\star}(k)$, can be viewd in Figure 6 (top right).

3.2 Alternating Direction Method of Multipliers (ADMM)

There are two natural approaches to solve the optimization problem (3). The first one is to calculate a centralized solution by optimizing the overall system at once. Here, a large rate of communication within the network is needed and the CE has to know all data of every prosumer, including possible charging rates, battery capacity, and current SOC, in detail. The second one renounces the communication by optimizing the behaviour of each prosumer separately. In this case the overall optimum cannot be guaranteed, see, e.g. [2]. A remedy is given by distributed optimization algortihms like dual decomposition or ADMM, see, e.g. [30] for further details.

ADMM is an algorithm used to solve large-scale optimization problems of the form

[TABLE]

in a distributed fashion. Here, the extended functions $f_{i}:D\to\mathbb{R}\cup\{\pm\infty\}$, $i\in[1:\mathcal{I}]$, $g:D\to\mathbb{R}\cup\{\pm\infty\}$, and $h:D\to\mathbb{R}\cup\{\pm\infty\}$, defined on the closed and convex domain $D$ are supposed to be proper, i.e. the effective domain $\Set{x\in D}{f(x)\in\mathbb{R}}$ is nonempty and there does not exist $x\in D$ such that $f(x)=-\infty$ holds (the definition analogously holds for $g$ and $h$).

Remark 3**.**

The functions $f_{i}$, $i\in[1:\mathcal{I}]$, may be used to incorporate individual goal of the subsystems in the optimization, see [17] for a detailed explanation.

Next, we outline the ADMM algorithm used in the following.

The great advantages of distributed optimization algorithms such as ADMM are scalability and plug-and-play capability. Since Step 1 in Algorithm 2 can be parallelized, the number $\mathcal{I}$ of subsystems does not affect the performance. Moreover, the specific data of the individual systems mentioned above is not known to the CE, which ensures plug-and-play capability.

3.3 ADMM yields the Solution of the Scalarized Problem

To apply ADMM in our framework to find efficient points of (2.4), we first scalarize the problem in order to obtain (3). As a second preliminary step, we introduce the auxiliary variables $\mathbf{a}_{i}=(a_{i}(k),a_{i}(k+1),\ldots,a_{i}(k+N-1))^{\top}\in\mathbb{R}^{N}$, $i\in[1:\mathcal{I}]$, and use the notation $\mathbf{a}=(\mathbf{a}_{1}^{\top}\ldots\mathbf{a}_{\mathcal{I}}^{\top})^{\top}\in\mathbb{R}^{\mathcal{I}N}$ to rewrite (3) as

[TABLE]

In ADMM, each prosumer first minimizes the augmented Lagrangian $\mathcal{L}_{\rho}:\mathbb{R}^{\mathcal{I}N}\times\mathbb{R}^{\mathcal{I}N}\times\mathbb{R}_{\geq 0}^{2N}\times\mathbb{R}^{\mathcal{I}N}\to\mathbb{R}$

[TABLE]

w.r.t. its own power demand for some $\rho>0$. Then, the CE minimizes $\mathcal{L}_{\rho}$ w.r.t. the aggregated power demand and the auxiliary variable $\mathbf{s}$ based on the aggregated optimal power demand provided by the prosumers.

The following theorem recalls convergence properties of ADMM, see [17, Theorem 3.1 and Subsection 3.2].

Theorem 4** (Convergence of ADMM).**

Let the augmented Lagrangian be defined by (9). Then, if there exists a saddle point of the unaugmented Lagrangian $\mathcal{L}_{0}$, i.e. there exists some $(\mathbf{z}^{\star},\mathbf{a}^{\star},\mathbf{s}^{\star},\lambda^{\star})$ such that

[TABLE]

ADMM converges for all $\mathbf{z}^{0}\in\mathbb{R}^{\mathcal{I}N}$ and $\Pi^{0}\in\mathbb{R}^{N}$ in the following sense:

1. (i)

The sequence $(\mathbf{z}^{\ell}-\mathbf{a}^{\ell})_{\ell\in\mathbb{N}}$ converges to zero. 2. (ii)

The sequence $\left((\kappa g(\bar{\mathbf{a}}^{\ell})+(1-\kappa)h(\mathbf{s}^{\ell})\right)_{\ell\in\mathbb{N}}$ converges to the unique minimum of (3.3). 3. (iii)

The dual variable $\lambda^{\ell}$ converges to the optimal dual point $\lambda^{*}$.

Firstly, we show that ADMM indeed yields the desired values.

Corollary 5**.**

For all $\kappa\in[0,1]$ the optimal values of $g(\bar{\mathbf{a}})=\left\lVert\bar{\mathbf{a}}-\bar{\zeta}\right\rVert_{2}^{2}$ and $h(\mathbf{s})=\left\lVert\mathbf{s}\right\rVert_{2}^{2}$ in the optimization problem

[TABLE]

are equal to those in (3.3).

Proof.

Theorem 4 (i) and (iii) imply that the (optimal) penalty term

[TABLE]

is constant, from which the assertion follows. ∎

Next, we work out that the optimal value of (3.3) also uniquely determines the optimal values of the functions $g$ and $h$.

Proposition 6**.**

For each $\kappa\in[0,1]$ there exists an optimal solution $(\bar{\mathbf{a}}^{\star},\mathbf{s}^{\star})$ of Problem (3.3). For $\kappa\in(0,1)$ the solution is unique. Furthermore, for $\kappa=0$ or $\kappa=1$ the optimal value $h(\mathbf{s}^{\star})$ or $g(\bar{\mathbf{a}}^{\star})$ is unique, resp. In particular, $\mathbf{s}^{\star}$ is unique for $\kappa=0$ and $\bar{\mathbf{a}}^{\star}$ is unique for $\kappa=1$.

Proof.

First consider $\kappa\in(0,1)$ and the modified problem (3.3). We show that $f_{\kappa}$ is strictly convex. Since $f_{\kappa}$ is twice differentiable, it suffices to show that its Hessian $H_{f_{\kappa}}$ is positive definite, see, e.g. [31, Subsection 6.4]. It holds

[TABLE]

and therefore

[TABLE]

Furthermore, since $\kappa\in(0,1)$, the Hessian $H_{f_{\kappa}}$ is positive definite and thus $f_{\kappa}$ is strictly convex. Hence, there is a unique minimizer $(\bar{\mathbf{a}}^{\star},\mathbf{s}^{\star})$ for (3.3).

Now consider $\kappa\in\{0,1\}$. Since $f_{\kappa}$ is at least convex in this case, there exists a unique optimal value $\kappa g(\bar{\mathbf{a}}^{\star})$ or $(1-\kappa)h(\mathbf{s}^{\star})$ and hence, $g(\bar{\mathbf{a}}^{\star})$ or $h(\mathbf{s}^{\star})$ is unique, resp. Furthermore, strict convexity of $h$ ($g$) yields uniqueness of $\mathbf{s}^{\star}$ for $\kappa=0$ ($\bar{\mathbf{a}}^{\star}$ for $\kappa=1$). ∎

Due to Proposition 6 we can reformulate (3.3) as

[TABLE]

where the optimization is w.r.t. the set of feasible function value pairs.

A fundamental observation on existence and uniqueness of solutions of (3.3) can be directly deduced from Proposition 6.

Corollary 7**.**

For each $\kappa\in[0,1]$ Problem (3.3) has an optimal solution. For $\kappa\in(0,1)$ the solution is unique.

4 Pareto Frontier

This section is dedicated to the Pareto frontier of the optimization problem (2.4). Firstly, we characterize the Pareto frontier and provide a brief sensitivity analysis w.r.t. the initial SOC and the size of the tube constraints. Then, we quantify the trade-off between the two objectives using the concept of proper optimality.

4.1 Determination of the Pareto frontier

We study the Pareto frontier of (2.4) for the open-loop optimal control problem to be solved at an arbitrary but fixed time instant $k\in\mathbb{N}_{0}$ within the MPC scheme presented in Algorithm 1. We illustrate our results for $k=0$ emphasizing that the numerical findings can be directly transferred to other time instants.

Since we scalarize (2.4) using a weighted sum, we are able to use the following characterization of efficiency given as Propositions 3.9 and 3.10 in [32].

Proposition 8**.**

Consider MOP (6) and its scalarization

[TABLE]

with $\mu_{i}\geq 0$, $i\in[1:m]$.

1. (i)

If $\mu_{j}>0$ for some $j\in[1:m]$ and $x^{\star}$ is an optimal solution of (11), then $x^{\star}$ is weakly efficient for (6). 2. (ii)

If $x^{\star}$ is a unique optimal solution of (11), then $x^{\star}$ is efficient for (6). 3. (iii)

Let in addition $\mathcal{M}$ be a convex set and let $f_{i}$, $i\in[1:m]$, be convex functions. If $x^{\star}$ is weakly efficient for (6) then there exist $\mu_{i}\geq 0$, $i\in[1:m]$, such that $x^{\star}$ is an optimal solution of (11).

The main result of this section, a characterization of the Pareto frontier of (2.4), is stated in the following proposition. Note that the optimal vector $\mathbf{s}^{\star}$ and the optimal vector $\bar{\mathbf{z}}^{\star}$ are unique for $\kappa=0$ and for $\kappa=1$, resp. (see Proposition 6). To find efficient points for $\kappa=0$ and $\kappa=1$, we consider the auxiliary problems

[TABLE]

and

[TABLE]

resp.

Proposition 9**.**

Let us consider the multiobjective optimization problem (2.4). Then, the following statements hold true:

1. (i)

The set of weakly efficient points of (2.4) is compact and connected. 2. (ii)

Consider $\kappa\in(0,1)$ and the unique optimal solution $(g^{\star},h^{\star})$ of (3.3). Then, there exists a unique weakly efficient point $(\bar{\mathbf{z}}^{\star},\mathbf{s}^{\star})$ of (2.4) with $g(\bar{\mathbf{z}}^{\star})=g^{\star}$ and $h(\mathbf{s}^{\star})=h^{\star}$. Moreover, $(\bar{\mathbf{z}}^{\star},\mathbf{s}^{\star})$ is efficient. 3. (iii)

Consider $\kappa\in\{0,1\}$ and an optimal solution $(g^{\star},h^{\star})$ of (3.3). Then, there exists a unique weakly efficient point $(\bar{\mathbf{z}}^{\star},\mathbf{s}^{\star})$ of (2.4) with $g(\bar{\mathbf{z}}^{\star})=g^{\star}$ and $h(\mathbf{s}^{\star})=h^{\star}$. If $g(\bar{\mathbf{z}}^{\star})$ or $h(\mathbf{s}^{\star})$ is the optimal value of Problem (12) for $\kappa=0$ or Problem (13) for $\kappa=1$, then $(\bar{\mathbf{z}}^{\star},\mathbf{s}^{\star})$ is efficient. 4. (iv)

The function $P:[0,1]\to\mathbb{R}^{2}$, $\kappa\mapsto(g(\bar{\mathbf{z}}^{\star}),h(\mathbf{s}^{\star}))$, where

[TABLE]

is well-defined and the component $P_{1}$ strictly decreases whereas $P_{2}$ strictly increases. Furthermore, $P$ is injective and the image of $P$ coincides with the Pareto frontier of (2.4).

Proof.

We show the single claims of the proposition consecutively.

1. (i)

As a subset of the compact feasible set of (2.4) the set of weakly efficient points is bounded. Moreover, Theorem 4.6 in [33] provides closedness and connectedness of the set of weakly efficient points of (2.4) and hence, compactness, which shows the assertion. 2. (ii)

Due to construction of (3.3) Statement (ii) follows directly from Proposition 8 (ii). 3. (iii)

Similar to the second claim the first part of Statement (iii) can be deduced from Proposition 8 (i). To prove the second part assume without loss of generality $\kappa=0$. Due to Proposition 6 there exists a minimizer $(\bar{\mathbf{a}},\mathbf{s}^{\star})$ of (3.3) and hence a minimizer $(\bar{\mathbf{z}},\mathbf{s}^{\star})$ of (3) where $\mathbf{s}^{\star}$ is unique. Moreover, since $g$ is strictly convex, the optimal solution $\bar{\mathbf{z}}^{\star}$ of Problem (12) is unique. By construction, $(\bar{\mathbf{z}}^{\star},\mathbf{s}^{\star})$ is efficient for (2.4). 4. (iv)

Well-definedness of $P$ is directly inherited from $g$ and $h$. The monotonicity properties can be seen as follows. For given $\kappa\in[0,1)$, choose $\varepsilon\in(0,1-\kappa]$. Then, we observe that $P_{2}(\kappa+\varepsilon)<P_{2}(\kappa)$ implies $P_{1}(\kappa+\varepsilon)>P_{1}(\kappa)$ and vice versa, i.e. if one component of the scalarized objective function strictly decreases, the other component has to strictly increase. Otherwise, this would directly yield a contradiction to the optimality of $P_{2}(\kappa)$ and $P_{1}(\kappa)$.

We prove the assertion by contradiction. To this end, let us assume that $P_{1}(\kappa)<P_{1}(\kappa+\varepsilon)$. Due to our preliminary considerations, this implies $P_{2}(\kappa+\varepsilon)<P_{2}(\kappa)$. Next, optimality yields the inequalities

[TABLE]

Due to our assumption and its implication, which show that the differences in brackets are strictly positive, we get

[TABLE]

a contradiction to $\varepsilon>0$.

Note that these considerations also imply the assertion for $\kappa=1$ since the case $P_{2}(1)<P_{2}(\kappa)$ for some $\kappa\in[0,1)$ is equivalent to $P_{2}(\kappa+\varepsilon)<P_{2}(\kappa)$, where $\varepsilon=1-\kappa$. This yields $P_{1}(\kappa+\varepsilon)>P_{1}(\kappa)$ and hence, $P_{1}(1)>P_{1}(\kappa)$, a contradiction to the optimality of $P_{1}(\kappa)$.

Since the components $P_{1}$ and $P_{2}$ are strictly monotone, $P$ is injective. Furthermore, Proposition 8 (iii) yields that each efficient point $(\bar{\mathbf{z}}^{\star},\mathbf{s}^{\star})$ of (2.4) can be found by varying $\kappa\in[0,1]$ and solving (3). Hence, by construction,

[TABLE]

which completes the proof.

∎

To determine the Pareto frontier numerically we run Algorithm 2 for different $\kappa\in[0,1]$. Here, we set \underaccent{\bar}{\mathbf{c}}_{i} and $\bar{\mathbf{c}}_{i}$, $i\in[1:\mathcal{I}]$, according to the values depicted in Figure 6 (top right). Let $\mathbf{z}_{i}(k)\in\mathbb{D}_{i}$, $i\in[1:\mathcal{I}]$, and $\Pi^{0}\in\mathbb{R}^{N}$ be given. For given $\kappa\in[0,1]$, the unique optimal values $P_{1}(\kappa)$ and $P_{2}(\kappa)$ are denoted by $g_{\kappa}$ and $h_{\kappa}$, resp. The results can be found in Figure 5. Note that for $\kappa=0.6$ the maximum of $f_{\kappa}$ is obtained.

Next, we display the values in the table in Figure 5 in a $(g-h)$-plane to get an approximation of the Pareto frontier, which is visualized in Figure 6 (top left). Note that for $\kappa\to 0$ and $\kappa\to 1$ the arcs asymptotically go to the minimum value of the optimization problem $\min_{(\bar{\mathbf{a}},\mathbf{s})\in\mathbb{S}}\,h(\mathbf{s})$ and $\min_{(\bar{\mathbf{a}},\mathbf{s})\in\mathbb{S}}\,g(\bar{\mathbf{a}})$, resp.

The convergence of the Pareto frontier in dependence of number of iteration steps within the ADMM-algorithm 2 is depicted in Figure 6 (bottom left). Here, we consider the functions $g:\mathbb{R}^{N}\to\mathbb{R}$ and $h:\mathbb{R}^{2N}\to\mathbb{R}$ as defined above and in addition we introduce functions $\hat{g},\hat{h},\hat{f}_{\kappa}:\mathbb{R}^{N}\to\mathbb{R}$ with

[TABLE]

Due to Theorem 4 the sequence $(\bar{\mathbf{z}}^{\ell},\mathbf{a}^{\ell},\mathbf{s}^{\ell})_{\ell\in\mathbb{N}_{0}}$ converges to the optimal solution $(\bar{\mathbf{z}}^{\star},\mathbf{a}^{\star},\mathbf{s}^{\star})$ of (3.3) fulfilling $g(\bar{\mathbf{a}}^{\star})=\hat{g}(\bar{\mathbf{z}}^{\star})$ and $h(\mathbf{s}^{\star})=\hat{h}(\bar{\mathbf{z}}^{\star})$. Regarding Figure 6 (bottom left), however, it can be observed that the convergence rate of $\bar{\mathbf{z}}^{\ell}\to\bar{\mathbf{z}}^{\star}$ is much higher than the convergence rate of $(\bar{\mathbf{a}}^{\ell},\mathbf{s}^{\ell})\to(\bar{\mathbf{a}}^{\star},\mathbf{s}^{\star})$. Note that the points generated by $g$ and $h$ (before convergence) are not feasible while those generated by $\hat{g}$ and $\hat{h}$ are not optimal.

The remaining graphic of Figure 6 (right) shows the open-loop perfomance (top) and the behaviour of the residual $|f_{\kappa}-\hat{f}_{\kappa}|$ depending on both the weighting parameter $\kappa\in[0,1]$ and the number of maximal iterations (bottom).

Remark 10**.**

Note that, due to Proposition 9 (v) the Pareto frontier as depicted in Figure 6 (top left) can be regarded as a the graph of the function $H:P_{1}([0,1])\to P_{2}([0,1])$, $H(P_{1}(\kappa))=P_{2}(\kappa)$, where $P$ is defined as in Proposition 9.

4.2 Sensitivity Analysis

In Figure 6 we set the average initial SOC $\bar{x}(0)=1.0$ and the tube size \bar{\mathbf{c}}-\underaccent{\bar}{\mathbf{c}}=0.2. Next we explore how the Pareto frontier changes if we modify these parameters.

In Figure 7 (left) different Pareto frontiers depending on the tube size are plotted. Increasing the size of the tube equals a relaxation of the subproblem $\min h(\mathbf{s})$. The larger the difference $\bar{\mathbf{c}}-\underline{\mathbf{c}}$, the easier it gets to keep $h_{\kappa}$ small. Small tube sizes do not affect the performance of $g$ drastically. If the tube becomes big enough to keep $h(\mathbf{s})$ small without effort, however, the focus of the optimization shifts to $g(\bar{\mathbf{a}})$.

The impact of the initial SOC $\bar{x}(0)$ on the Pareto frontier is visualized in Figure 7 (right). At first sight, the performance w.r.t objective $g$ seems to improve with increasing $\bar{x}(0)$ while the performance w.r.t. $h$ deteriorates and vice versa. To interpret this behaviour one has to take the time series $\bar{w}(n)$, $n\in[k:k+N-1]$, into acount, see Figure 6 (top right, dotted, black line). If the batteries are completely discharged at the beginning of the time interval there is no possibility to reduce the energy demand by discharging the batteries. Hence, the tube constraints are violated due to $\bar{\mathbf{w}}(n)>\bar{c}(n)$ for small $n$. To compensate $\sum_{n=k}^{k+15}\bar{\zeta}(n)>\sum_{n=k}^{k+15}\bar{\mathbf{w}}(n)$, on the other hand, one needs to charge the batteries. Hence, the higher the SOC at the beginning, the harder it gets to trace the desired trajectory.

4.3 Proper Optimality

Considering an efficient point, improving the performance w.r.t. one objective is only possible at the expense of the performance w.r.t. the other. We investigate the trade-off between (potentially) conflicting objectives using the concept of proper optimality. Note that there are different definitions given by several authors [34, 35, 32]. In this paper, however, we focus on proper optimality in the sense of Geoffrion, which expands efficiency by introducing an upper bound on the trade-off between two objectives.

Definition 11** (Proper efficiency).**

A point $x^{\star}\in\mathcal{M}$ is called properly efficient (in the sense of Geoffrion) of the MOP (6) if $x^{\star}$ is efficient of (6) and there exists some $L>0$ such that for all $i\in\{1,\ldots,m\}$ and $x\in\mathcal{M}$ with $f_{i}(x)<f_{i}(x^{\star})$ there exists some index $j\in\{1,\ldots,m\}$ with $f_{j}(x^{\star})<f_{j}(x)$ and

[TABLE]

Remark 12**.**

In Figure 8 (top right) we consider the optimal solution $(g_{\kappa_{0}},h_{\kappa_{0}})$ of (3.3) corresponding to a given $\kappa_{0}\in[0,1]$. If we want to improve the perfomance w.r.t. objective $g$ by the amount of $g_{\kappa_{0}}-g_{\kappa_{1}}$, the performance w.r.t. objective $h$ worsens by the amount of $h_{\kappa_{1}}-h_{\kappa_{0}}$. Using function $H$ as defined in Remark 10 and the relation $g_{\kappa}=P_{1}(\kappa)$, $\kappa\in[0,1]$, Inequality (15) can be written as

[TABLE]

For the particular case $g_{\kappa_{1}}\to g_{\kappa_{0}}$ this reads as

[TABLE]

where the left hand side is the negative of derivative of $H$ w.r.t. $g_{\kappa}$ at $g_{\kappa_{0}}$ if existent. The latter indicates the cost of an infinitesimal improvement w.r.t. the performance of $g$. Hence, proper optimality provides an upper bound to this cost which is an essential information for a decision maker.

The following theorem provides a characterization of the properly efficient points of (6) using a weighted sum scalarization as in Proposition 8.

Theorem 13** ([32], Theorem 3.15).**

Consider MOP (6) with a convex feasible set $\mathcal{M}$ and convex component functions $f_{i}$, $i\in[1:m]$. Then the following holds true: A point $x^{\star}\in\mathcal{M}$ is properly efficient for (6) if and only if $x^{\star}$ is an optimal solution of (11) with positive weighting parameters $\mu_{i}>0$, $i\in[1:m]$.

This result can be directly applied to (2.4).

Corollary 14**.**

All efficient points of (2.4) except for the extremal points $(\kappa\in\{0,1\})$ are properly efficient.

We conclude this section with an example of calculating the bound $L$ for $\kappa=0.2$ numerically.

Example 15**.**

Consider $\kappa_{0}=0.2$. We provide an intuitive way to determine an upper bound $L$ at the corresponding optimal value pair $(g_{\kappa_{0}},h_{\kappa_{0}})=(1.452,0.442)$. Since $g_{\kappa}$ and $h_{\kappa}$ result from an optimization routine for which there is no explicit formula, we cannot compute $L$ explicitly.

Instead, we calculate

[TABLE]

for all $\kappa_{1}$ in the table in Figure 5 (left). Additionally we compute $L_{\kappa_{0}}(0.01)$ and $L_{\kappa_{0}}(0.99)$ to illustrate the behaviour near the extrema. Then, we choose $\max\{L_{\kappa_{0}}(\kappa_{1})|\kappa_{1}\in K\}$, where $K=0.05\cdot[0:20]\cup\{0.01,0.99\}$, as an approximation of the upper bound $L_{\kappa_{0}}$ in (15). The results can be found in the table in Figure 8.

One can observe that $L_{0.2}\approx 4.74$ holds, i.e., to gain one quantity in $g$ direction one has to spend at maximum 4.74 quantities in $h$ direction and vice versa. In our case, i.e. $\kappa_{0}=0.2$, the bound on the second relation is active while the first one can replaced by $1/L_{0.2}=0.211$, i.e., to gain one quantity in $g$ direction one has to spend at maximum 0.211 units in $h$ direction.

In Figure 8 (bottom right) the evolution of the bound $L_{\kappa}$ is illustrated depending for $\kappa\in[0.05,0.99]$.

5 Conclusions and outlook

In this paper, we have analysed the Pareto frontier of a multiobjective optimization problem, in which a trade-off between peak shaving and providing flexibility has to be made. To this end, we first considered the corresponding scalarization, linked it to ADMM in order to (numerically) compute its solution for a given scalarization parameter, and rigorously showed that by doing so we get the whole Pareto frontier. Moreover, by using the concept of proper optimality introduced by Geoffrion, we further investigated the Pareto frontier, which allowed to quantify the trade-off between the conflicting objectives.

There are several ways to extend our model: One could easily consider the impact of controllable loads as introduced in [36] or varying lower/upper bounds on (dis)charging rates. Another possible modification of our optimization problem is to additionally penalize the use of the batteries, e.g. in form of $\left\lVert u_{i}\right\rVert$. The presented approach, especially ADMM still works, as long as the sets $\mathbb{D}_{i}$ of feasible solutions $\mathbf{z}_{i}$, $i\in[1:\mathcal{I}]$, stay polyhedral.

Throughout this paper, we assumed the net consumption to be given. Finding a suitable prediction method is an interesting topic for future research. For the goal of load shaping, the robustness of MPC schemes w.r.t. inaccurate forecasts was numerically investigated in [2]. Moreover, there are interesting approaches to estimate the energy generation within the prediction window more thoroughly using artificial neural networks as pointed out in [37, 38]. A good starting point would be [39] and the references therein.