Model Predictive Frequency Control
Employing Stability Constraints
Abstract
In this article a model predictive control (MPC) based frequency control scheme for energy storage units was derived, focusing on the incorporation of stability constraints based on Lyapunov theory and the concept of passivity.
The proposed control schemes, guaranteeing closedloop stability, are applied on a onearea and twoarea power system. For the twoarea power system, a coordinated (centralized) control and an uncoordinated (decentralized) control approach is conducted.
The stability properties of the different MPC setups were derived, implemented and simulated. Furthermore the corresponding control performance was analyzed.
I Introduction
Ia Motivation
Traditionally, power system operation has essentially been based on the assumption, that electricity is reliably and steadily produced by large power plants, which are fully dispatchable, i.e. controllable and have a high frequency inertia.
However, a strong trend towards generation of electrical energy by renewable energy sources (RES), i.e. PV units and wind turbines with no or decoupled rotating masses, exists^{1}^{1}1On 16th of June 2013, the share of energy produced by wind and solar power reached a level of above 60 % for the first time in Germany. It resulted in negative electricity prices in France and Germany [1].. How power systems could be adapted in order to accommodate for increasing shares of uncontrolled fluctuating RES as well as power market activities is a highly relevant and still open research question.
Options to deal with these challenges might be the expansion of transmission capacities, the extensive integration of storage capacities as well as a better exploitation of the inherent flexibility within the power system.
Model predictive control (MPC) as an optimal control scheme for regulating grid frequency receives rising attention due to rapidly growing shares of variable RES and thereby arising challenges for power system operation.
The choice of MPC as a control approach is especially motivated by its ability to incorporate operational constraints of power systems, which cannot be handled by conventional P/PIcontrollers. It enables the provision of frequency control using a generic power system storage unit, e.g. a battery, with given operational constraints, such as the power ramp rate, power rating and storage capacity.
Additionally, the recently growing interest in using MPC for control purposes of power systems emerges due to the availability of faster and cheaper computing resources, as a significant computational effort comes with the use of such a receding horizon control scheme.
IB Goal and Content
The objective of this article is to derive an MPC based frequency control scheme for energy storage units. Specifically, the employment of stability constraints will be studied. The theoretical focus lies on the derivation of passivity as a concept for guaranteeing stability. Furthermore stability guaranties derived from LyapunovTheory, so called Control Lyapunov Functions (CLF) are incorporated into the optimal control setup.
Eventually, the goal is to ensure a stable frequency regulation of a onearea and a twoarea power system, which is disturbed by a fault signal.
The remainder of the article is organized as follows: Section II elaborates on the dynamics of grid frequency. Subsequently in section III, first, a theoretical introduction to stability is stated, followed by derivation, simulation and validation of conventional control.
Subsequently in section IV, the derivation of the model predictive frequency control (MPFC) problem is given. In section V different approaches to ensure stability of the controller are derived. This is followed in section VI by a study case. Finally, the outcome of the article is presented in section VII.
Ii Power System Modeling and Analysis
Iia Grid Frequency
Power systems are dynamic systems with a high degree of complexity. Corresponding processes within the power system therefore happen on several timescales (milliseconds to yearly seasons), with a large spatial distribution (hundreds to thousands of kilometers) and numerous grid hierarchy levels (voltage levels).
Dominant for the effect of frequency stability, and hence stable operation, is active power balance, meaning that the total power infeed minus the total consumption is zero. The nominal frequency is set to Hz as in the ENTSOE Continental Europe system (formerly UCTE).
Generally deviations from the nominal frequency arise from imbalances between instantaneous electric power consumption and production. If there is a higher consumption than production in a time instant, this results in a decelerating effect on the synchronous machines. On the contrary, higher production than consumption of electric power, leads to an accelerating effect on the synchronous machines [2].
Small local disturbances can evolve into consequences influencing the whole power system, which could lead to cascading faults and blackouts in the worst case.
The key system state to observe is the grid frequency , concerning this (dominant) aspect of power grid stability. As the rotational speed of the synchronous generators is inherently coupled to grid frequency (), frequency deviations should be kept as small as possible, since this leads to damaging vibrations in synchronous machines.
Maintaining the grid frequency within an acceptable range is therefore required for stable operation of the power system. Small variations occur spontaneously by small load or generation deviations and usually do not have critical consequences in normal operation.
Large frequency variations, which might be caused by errors in demand forecasts or RES forecasts, the spontaneous loss of load or generation units, however, could lead to situations, where synchronization between the generating unit is lost. This happens when the angle differences get too large, which would result in (possibly cascading) disconnection of the machines [3].
IiB Aggregated Swing Equation
The frequency dynamics of a power system can be described by the aggregated swing equation (ASE)
(1) 
with being the deviation of load, being the mechanical (or inert) power deviation, and being the total inertia constant^{2}^{2}2 is typically valued in the range 2…10 s, cf. [4]..
Please note, that losses are not taken into account here. Regarding the stability analysis, which is done in the sequel, this is a very conservative assumption.
The following assumptions are made

The synchronous machine is modeled as a constant electromagnetic field behind the transient reactance. The angle of the electromagnetic field is assumed to coincide with the rotor angle.

Resistances in lines, transformers and synchronous machines are neglected.

Voltages and currents are assumed to be perfectly symmetrical, i.e. pure positive sequence.

The angular velocity is close to nominal.

Static models for lines are used.
Despite its simplicity, a number of important conclusions concerning the angular stability in large systems can be drawn from this analysis [5].
IiC Frequency Dependency of the Loads
In real power systems, a frequency dependency of the aggregated system load is clearly observable. This has a stabilizing (”selfregulating”) effect on the system frequency .
It is assumed, that the damping power of the system can be written as
Including the frequencydependency of the loads into equation 1 leads to
(2) 
Rotating mass loads are neglected. They play a decreasing role in future power systems and are not of big importance for (conservative) stability analyses.
IiD Mathematical Modeling
For the system dynamics the following basic model is assumed:
(3) 
where with referring to the deviation of the normal frequency of the grid and to the state of charge of a connected battery.
The corresponding linearization in as a state space model is given by:
(5) 
Next, the continuous linear system is discretized
(6) 
To be able to account for energy constraints of the control input later (due to energy storage restrictions etc.), the additional power is integrated as following, and is used as another state within the state space model:
Including the possibility of losses, respectively selfdischarging of the battery with capacity leads to
(7) 
with . Combining equation 7 and the aggregated swing equation 2 and transforming it into (continuous) state space form yields:
(8) 
with:
refers to the deviation of the expected power, including the fault power, therefore, it directly influences the frequency deviation in the grid.
however refers to the electrical power, which is injected into the grid by the controller, hence directly influencing the state of charge of the battery [6], [7].
IiE TwoMachine Model
If a power system is highly meshed, it can be seen as all units being connected to one single bus, i.e. acting as one (aggregated) swing equation. However, in practice, a large interconnected system is divided into several ”control zones”. To understand the interaction between two such zones, in this section a twoarea model as in figure 1, with a tieline in between is derived.
The power exchange between the two zones is
where equals the (equivalent) reactance of the tie line between zone 1 and zone 2. This can also be written as:
with being the maximum power transmission
The (aggregated) swing equation divides into a system of equations, comprising the interaction between the two areas:
(9) 
After linearizing the model around the set point , where , the state space matrices change to
and
with
the corresponding states and inputs
IiF Simulation of Power System Dynamics
For simulation purposes, power system parameters as in table I were assumed. For load damping, the minimum measurement in [8] was assumed, to carry out very conservative stability analysis.
Unit  Size  Description 

6 s  system inertia  
self regulating effect  
load damping  
frequency droop, for primary freq. control  
240 s  integration param. for conv. freq. control  
0.17  proportional param. for conv. freq. control  
= 20550  bias factor, for conv. control 
IiF1 Fault Signal
In all simulations, as long as not stated otherwise, the disturbance signal given in the power deviation graph in figure 2 is introduced in the power system. This signal was taken as it provokes instabilities for (conservative) stability tests due to enormous sudden disturbances, e.g. cascading disconnection of load or generation units. Applying the fault to an uncontrolled, system comprising two control areas shows the frequency response in figure 3.
Iii Conventional Frequency Control
Iiia Power System Stability
A dynamic phenomenon in a power system is initiated by a disturbance in the system, e.g. a a line or a generator trips. While smaller disturbances usually result in small transients in the system, which diminish quickly because of damping, larger disturbances could have a significant impact and excite larger oscillations. Hence, stability is associated with decaying oscillations and that the operation of the power system can continue without major impact for any power consumer [5].
In this article the focus lies on frequency stability, corresponding to the frequency dynamics, which were analyzed in section II. Unlike voltage stability, which only acts locally, frequency stability plays a global role in a power system. It refers to the difference between the total power fed into the power system and the total power consumed by the loads, including the losses. If the resulting imbalance is comparatively small, the generators participating in frequency control will regulate the active power input from their prime movers, and bring back the frequency deviation to acceptable values. But if the imbalance is too large, the frequency deviation gets significantly big and might cause serious consequences.In real power systems, instabilities result both from active and reactive power imbalances, so the assumptions, which are made in this article, are of course not valid in every case of power system instabilities. However, in many cases it is possible to identify active power imbalances and resulting frequency instability as the dominating processes at the start of power system instabilities [5].
IiiB Conventional Frequency Control
The nominal grid frequency is defined to be Hz in the grid zones of the European Network of Transmission System Operators for Electricity (ENTSOE). To maintain this frequency, for secure operation, electricity demand and supply have to be in balance at every point in time. Small local disturbances in a power system can lead to consequences influencing the whole system, in the worst case ending in a blackout.
In order to secure a stable operation, a frequency regulation mechanism is implemented. In traditional frequency control, there are three categories (cf. figure 4): Primary frequency control provides power output proportional to the deviation of the system frequency, which stabilizes the system without restoring the reference frequency. Its time horizon lies within a few seconds after the occurrence of the deviation. Furthermore secondary frequency control, taking over after approx. 30 s, has an integral control part, which restores the reference frequency. The time constant of secondary frequency control is chosen to be significantly larger than the time constant for primary control in order not to interfere with it as well as avoid ”wear and tear” for the units. Tertiary control, finally, is operated manually and adjusts power generation setpoints about 15 min after severe faults. Generator rescheduling is dispatched through intraday auctions, according to the estimated permanent fault magnitude, aiming to relieve tertiary control by cheaper sources [3].
One of the factors used in this control scheme is the frequency bias factor. Load frequency control is based on the noninteraction principle, hence, the disturbance power balance between all neighboring areas is restored [9].
IiiC Simulation of Conventional Frequency Control
One Area
Two Areas
In this setup, two equally sized areas are interconnected via a tieline, e.g. with a transmission capacity of p.u. Each controlled area is equipped with an individual conventional primary and secondary controller.
IiiD Motivation for new Methods
In conventional systems, primary and secondary frequency control, referred to as automatic generation control (AGC), usually is implemented with satisfactory results.
The main conventional AGC properties are (summarized from [11]):

”The control system is naturally areawise decentralized.” Except in the case of a lack of area generation capacity, the frequency deviation diminish and all the interarea power transmissions return to their preassigned values, while each area needs only data available within the area itself.

”The control system is robust.” Robustness is meant in a ”rustic, simple and not fragile way”:

Only regional data is needed.

Only little information is needed: frequency, the sum of exported powers and the electric generated powers, with possibly low sampling rates and large transmission delays.

It is practically insensitive to the exact physical system parameter values as well as control system failures.


”The control system is sufficiently stable in practice,” if it is carefully designed and tuned. Both static and dynamic aspects have been studied in the past.

”Adaptive control is easy to implement.” This is due to the system’s simplicity.

”Control may not be sufficiently smooth.” No attention is paid to unit response rate limits in the controller. In practice, regulating margins are used, but rate constraints can not be guaranteed to be met.

”The interface between load frequency control^{3}^{3}3I.e. primary and secondary control. and economic dispatch^{4}^{4}4I.e. tertiary control and generation rescheduling. is not satisfactory.” Due to different time scales, the orders of the different control levels might contradict each other, resulting in undesired oscillations as well as useless wear of the generation units.
Besides the final two properties, the given points essentially explain the success and long life of the conventional implementations.
But, many proposals exist for modern control systems, most of them applying optimal control theory or optimal power flow techniques. The relevant capabilities and desirable properties of these modern proposals are seen as high (summarized from [11]):

It could take system constraints into account, e.g. transmission security and power rate limits.

It might improve economy, by an accurate application of economic dispatch and the reduction of regulating margins.

It still must have good transients, with large stability margins and smooth control, but especially addressing the LFCED interface problem.

It must be kept robust, i.e. not fragile, by avoiding data complexity, and, probably most important, being areawise decentralized.
With all potentials lying in modern control techniques involving optimal control and optimal power flow: one should never forget the engineer’s perspective, particularly regarding the objectives, the constraints and, especially, the robustness of the control system [11].
Iv Model Predictive Control
Iva Main Concept of Model Predictive Control
The design of a stabilizing feedback with a performance criterion being minimized, while satisfying constraints on the control input and the system states at the same time, is desired for many control problems. A solution approach is to repeatedly solve an open loop optimal control problem. The first sample of the resulting openloop optimal trajectory, is then applied to the real system. Afterwards the whole process is repeated for the next time step. Control methods like these are referred to as either model predictive control (MPC), moving horizon or receding horizon control (RHC) [12].
MPC is formulated as a repeated online solution of a finite horizon openloop optimal control problem. When measuring the initial system state , a control input over the prediction horizon can be determined, by solving the following optimal control task:
s.t.  (10)  
with the control input and the system state . is an integer counting in time: , refers to the cost at time , to the state at time , to the corresponding, cost minimizing control input at time while is the (allowed) solution space for the control inputs and the respective (allowed) solution space for the system states .
IvB Model Predictive Frequency Control
The model predictive frequency controller uses the directly available system output as control signal. The control input (determined by the controller) is the regulation power that is supplied to the power system (e.g. by a battery).
The system, implemented within the MPC controller, based on the model in equation 8 derived in section II, is the discretization of:
(11) 
with
The first line in equation 11 refers to the swing equation, while the second line captures the state of charge by integrating and accounts for losses. The MPC cost function weights^{5}^{5}5The ramp rate is not penalized, as this not of substantial interest for energy storages like batteries, that can ramp power output very fast. are defined as , , with steps and a sampling time of ms:
s.t.  (12)  
V Employing Stability Constraints
Va Control Lyapunov Function
Derived from Lyapunov theory, another approach for guaranteeing stability can be found. The new approach is based on a finite horizon optimal control problem with terminal cost. By inserting a socalled Control Lyapunov Function (CLF) as a terminal cost term
into a generalized MPC setup, stability is achieved.
For a nonlinear system
where , and
a CLF is a proper, positive definite function such that:
If it is possible to ensure a negative derivative at every point by an appropriate choice of , the system can be stabilized with . It can be shown, that such a CLF exists, when a globally asymptotically stabilizing control law exists (smooth everyhwere except at , cf. [13]). However, there do not exist systematic approaches to find suitable CLFs for given nonlinear systems [14]. The terminal region is not restricted .
VA1 CLF based MPC
To implement a Control Lyapunov Function as a stability guarantee into an MPC setup, the corresponding terminal cost has to be derived. For that, a quadratic penalization of the form is used: , with being the solution of the Lyapunov equation:
(13) 
The solution is symmetric and positive definite, as is symmetric, positive definite and has all its eigenvalues inside the unit disk (time discrete case). This approach assumes no control action after the end of the horizon, so that
This only makes sense if the system is asymptotically stable, or no solution will exist.
VB Passivity
As with the Zero Terminal State constraint and the Control Lyapunov Functions described in the preceding part, closedloop stability can also be guaranteed by imposing passivity for a model predictive control problem. In [15] a nonlinear model predictive control scheme was proposed, which is based on the optimal control, nonlinear model predictive control and Control Lyapunov Function. Analogously in [16] a nonlinear model predictive control scheme was developed, which makes use of the relationship between passivity and optimal control, as well as the relationship between nonlinear model predictive control and optimal control.
In the sequel, an introduction into the concept of passivity is given. Passivity is tied to optimal control by an input affine system, which is optimal if and only if it satisfies a passivity property with respect to the optimal feedback [17].
Based on this, it is shown in [16], that passivity and (nonlinear) model predictive control can be merged together, such that the individual advantages of each concept can be maintained, i.e. feasibility and closedloop stability due to passivity as well as good performance due to online optimization in nonlinear model predictive control.
Following [16], consider the following affine nonlinear system
(14) 
where is the state, is the input and the output. Local Lipschitz continuity as well as being an equilibrium point is assumed. Then, system 14 is said to be passive, if there exists a positive semidefinite storage function , such that the following is satisfied for all :
(15) 
where is a solution of 14. If is differentiable as a function of time, then the this relation can be reduced to
(16) 
In case, system 14 has a welldefined normal form, a further characterization of passive systems is also possible in terms of relative degree and minimumphase property. Specifically, the relative degree needs to be , i.e. , and the system must be weakly minimumphase.
VB1 Passivity and Stability
To achieve asymptotic stability of system 14, one can make use of the relationship between Lyapunov stability and the concept of passivity. It can be established by using the storage function as Lyapunov function. As passivity only requires to be positive semidefinite, the equlibrium point might be unstable, even if passivity is ensured. This would be the case, if there exists an unobservable part of the system that is unstable. Hence, this unobservable part needs to be asymptotic stable for system 14, i.e. for its solution holds:
(17) 
which is is the so called property of zerostate detectability.
In [18] a concise summary of the key points describing the relationship between passivity and Lyapunov stability is given:

The equilibrium point, , of the system 14 with zero input, , is stable, if the storage function is positive semidefinite and the system is passive and zerostate detectable.

The equilibrium point, , of the system 14 with zero input, , is stable, if the storage function is positive definite and the system is passive.

The equilibrium point, , of the system 14 with zero input, , is asymptotically stable, if the storage function is positive definite and the system is strictly passive.
In this case, a semidefinite storage function is considered. Assuming zerostate detectability, the system can be stabilized with the feedback . With equation 16, the following relation can be established:
(18) 
and [17]. Furthermore, if is radially unbounded, i.e. for , and all solutions are bounded, then the system 14 can be globally stabilized by feedback .
Please note here, that 18 can promptly be rewritten to an inequality, which easily adds to common MPC setups as a state constraint:
(19) 
Unfortunately, the derived passivity constraint is generally not convex, written for a timediscrete setup:
with and written for the case with one state and one control input :
where has the eigenvalues and , which means, that is not positive (semi)definite. Actually, by a similar argument, is bilinear and therefore never convex.
To obtain a convex QP optimization problem, which can be solved by standard solvers, the passivity constraint is implemented in a way, such that it only holds for the first (applied) sample of the receding horizon. With this, is fixed, a linear and hence convex constraint is obtained:
(20) 
Passivity and Optimality
In addition to the statements given in the sections before, the relationship between optimal control and passivity can be established by using the value function as a storage function and the otpimal feedback as an output of the system 14. Taking the optimal feedback, which stabilizes the considered system, it minimizes the performance index if and only if the system
(21) 
is zerostate detectable and output feedback passive with respect to
VB2 Passivity based MPC
The in the preceding parts of this article, both (nonlinear) model predictive control and passivity were presented. Here, the aim is to merge these two concepts, inspired by the relationships between optimal control and passivity and between optimal control and nonlinear model predictive control [15].
The goal of merging those concepts is to maintain their corresponding individual advantages, i.e. good control performance due to online optimization within the nonlinear model predictive scheme and guaranteed closedloop stability due to passivity, referred to as passivity based NMPC.
Suppose a general NMPC setup, with a passivity constraint added (last line):
(22)  
Where the last line directly follows from 19 as a passivity based constraint. In case the system 14 is passive and zerostate detectable, it can be stabilized with the feedback (cf. [15]). Therefore, the passivity based constraint is a stability constraint which guarantees closedloop stability. This holds even globally, if the storage function is radially unbounded and all solutions of the system are bounded.
Unlike to many other nonlinear model predictive control schemes, which achieve stability by enforcing a decrease of the Control Lyapunov Function (CLF), , along the trajectory of the solution, stability is achieved directly with the aforementioned state constraint.
On first sight, the proposed scheme seems to be rather restrictive, as it is only applicable to passive systems. However, for stabilizing purposes, no real physical output is needed. It is enough to have a fictitious output , yielding the following (fictitious), passive system:
(23) 
Note, that there always exists such a fictitious output , if a Control Lyapunov Function exists, since then, by definition is a (fictitious) passive output. Unfortunately, there is no standardized way to construct a passive output. But as passivity is a physically inherited concept, such a fictitious passive output can often be found.
Applied Passivity Constraint
In the following, the applied passivity constraint will be derived. As stability with respect to the frequency deviations and should be secured, a storage function of the form is taken:
(24) 
The goal is, to fulfill the passivity property
(25) 
with
(26) 
First, the constant is derived, which is a proportionality factor of the frequency , describing the kinetic energy stored in the rotating machines of area around the setpoint :
with
(27) 
Assuming two machines or areas , the storage function from equation 24 becomes:
(28) 
where
is true for all , since and
.
The passivity property (cf. equation 26) hence yields
(29) 
Combined with equation 27, and as well neglecting all (marked) negative terms, equation 29 can be transformed into
(30) 
Taking yields and together with the passivity property (cf. equation 25) finally results in
(31) 
as the corresponding passivity constraint for a twomachine system. Analogously for the onemachine system:
(32) 
In the here presented implementation, the passivity based constraint is enforced only during the first sampling interval and not for the whole length of the finite prediction horizon . This allows the same stability guarantees as enforcement of the constraint for the whole horizon does, while it generally exhibits better feasibility and less computational effort (cf. section VB).
VC Summary
The passivity based nonlinear model predictive control approach combines passivity with model predictive control by demanding the derivative of a storage function to fulfill (being positive semidefinite). At the same time their respective advantages are maintained: Feasibility and closedloop stability due to passivity and good performance due to online optimization by model predictive control.
A stability guaranteeing Control Lyapunov function ensures a storage function to satisfy (which actually needs to be positive definite). Where such a Lyapunov function cannot be found, the passivity based approach is especially appealing.
The corresponding basic principles are the relationships of passivity and nonlinear model predictive control to optimal control.
Vi Study Case and Results
Via General Setup
The application analyzed for the derived stability concepts is to guarantee frequency stability in a power system. In this article, the focus lies on the evaluation of the interaction between two areas, interconnected with a tieline. For this case the used model was derived in section IIE.
All examples are modeled as a quadratic program (QP) and in perunit system (p.u.) for scalability. The primary optimization objective is to maintain a frequency deviation , meaning, that the frequency should be close to Hz in every area. Internally, the frequency related states of the used model are normed to . For illustrative purposes, in all shown diagrams, frequency related states are transformed to frequencies, hence they are shown with the unit Hz.
For simulation and optimization purposes, the YALMIP environment was used, developed and maintained by Johan Löfberg. The software package used for mathematical optimization within the predictive control framework of YALMIP is CPLEX. All simulations were run on a computer with a quadcore CPU (Intel Core i7, 3.4 GHz), with 16 GB memory and 64bit Windows 7 as operating system.
ViB Developed Simulink Model
A model was developed within Simulink, which is part of a MATLAB package.
An overview of the devloped Simulink model is provided in figure 8.
ViB1 One Area
For carrying out experiments in a one area setup, the control input for the second plant was disconnected and a zero transmission power () was assumed. The model equations used are time discretized versions of the derived model in section IID. The following state vector and input was used:
The MPC cost function weights^{6}^{6}6Please note at this point, that internally the frequency state penalized is normed to . are defined as:
with a prediction horizon length of steps and a sampling time of s. Observe here that penalizing with has physical and/or economical reasons, but might lead to instabilities using the standard MPC setup. The following setup (derived in section IVB) is referred to as the standard MPC setup:
s.t.  (33)  
For ensuring wide feasibility of the optimization algorithm, the limit of the frequency deviation was set to a rather big range (higher deviations might cause serious system damages): The state of charge is limited to: Initially, both states are assumed to equal zero. The battery power is constrained with
A rate constraint is implemented as well. It is assumed, that a battery can change the power by 1 p.u. at each sample, hence
Setting up Passivity
Setting up the CLF
To setup the CLF based MPC, the terminal costs for the one area system are defined as the solution of the (discrete) Lyapunov equation^{7}^{7}7In Matlab: q_term = dlyap(A(1,1),Q(1,1))., which results in This is one dimensional, as it only refers to stability of the first state, the frequency deviation. Stability of the second state, the state of charge, is implied by constraining the integrated control action. Therefore, the corresponding terminal costs for the second state are set to zero:
ViB2 Two Areas (Uncoordinated Case)
The uncoordinated two area system consists basically of two equal one area systems (each one comprising the same MPC setup as in the one area case), where power transmission between the two systems is activated. Both systems are controlled individually with an own energy storage system, therefore, all parameters and constraints are set as in the one area system mentioned before.
The power transmission is implemented as derived in section IIE in the form of:
with a maximum power transmission of
Setting up Passivity
Passivity is set as a constraint for each system individually:
with and it follows applied for the first sample:
Setting up the CLF
The CLF setup is done as in the one area control case (cf. section VIB1), individually for each controller.
ViB3 Two Areas (Coordinated Case)
In the coordinated controlled version of the two area system, the power angle between the two systems is introduced as a new system state. The corresponding model equations as well as the new state space matrices reflecting the coupling between the systems were derived in section IIE. Still, it is assumed, that every area uses its own energy storage unit, hence the two area system consists of the state vector and input vector
The cost function weights hence change and were defined as the following for the coordinated case:
Besides that, all parameters and constraints are the same as in the uncoordinated two area system.
Setting up Passivity
The passivity constraint is adapted to account for both frequency related states and and the same time:
respectively:
Setting up the CLF
For the coordinated version, CLF stability for both frequency deviations is desirable at the same time, hence the solution of the (discrete) Lyapunov equation becomes the terminal cost as in the following matrix, accounting for both frequency stability related states:
The corresponding terminal costs for the other states as the angle and the states of charge are set to zero:
ViB4 Summary
To summarize, besides the standard MPC setup (as in equation 33) the applied stability approaches are implemented as following (generalized for both one and two area control).
Standard MPC Setup
s.t.  (34)  
Passivity based MPC Setup
As an additional constraint to system 34, the passivity constraint is incorporated:
CLF based MPC Setup
Here, the terminal costs are adjusted, i.e. compared to the standard setup in system 34, the cost function changes to
ViC Results
Two different setups were analyzed. First, the uncoordinated case, where a decentralized control scheme was applied. Second, the coordinated case, where a centralized controller has knowledge about all (potentially distributed) system states.
First, the angle differences between the two systems are presented for a fixed prediction horizon . Afterwards, the simulations were generalized for various different prediction horizons .
The applied fault signal is the one shown in figure 2 for every simulation carried out. Note, that the error is not zero mean and becomes constantly zero at s.
ViC1 Two Areas (Uncoordinated Case)
Employing frequency stability for two equally sized, interconnected areas of power systems with a maximum transmission capacity of p.u. was simulated.
The uncoordinated version was analyzed with respect to standard MPC control, passivity based MPC control and CLF based MPC control. Uncoordinated means here, that each area is controlled in an individual way and each controller has only local information, hence, a decentralized control topology. This implies a significant modelplantmismatch, as the possibility of power transmission is not incorporated within the MPC model.
The effect of stability can be seen especially illustrative in figure 10, which shows the angle difference for the two controlled zones with respect to the different control approaches. Large angle differences (as conventional control exhibits) as well as increasing oscillations (as in standard MPC) indicate frequency instability.
ViC2 Two Areas (Coordinated Case)
The power system itself is not altered. However, the area coupling is incorporated into the MPC setup in the coordinated approach. By that, the modelplant mismatch is significantly reduced compared to the uncoordinated case. Coordinated control refers to a centralized control topology, where the controller has full information about all states in all control areas. Compared with the uncoordinated case, the stability behavior of a standard MPC approach improves, yet still exhibiting undesired persisting frequency oscillations. Essentially, the stability behavior of the different approaches remain unchanged.
ViC3 Generalizing Two Area Control
The behavior of the states with respect to time was shown precedingly for a fixed prediction horizon . Now, a generalization of controlling two areas for various prediction horizons is made.
For facile comparisons, the performance measures of the presented MPC approaches with employed stability constraints for two uncoordinated areas are shown together with the performance measures for two coordinated areas in figures 12, 12, 14, 14, 16, and 16. In these figure, control performance measures for two controlled systems (one being disturbed) over different horizon lengths, for two coordinated (continuous lines) and two uncoordinated areas (dashed lines) are compared.
Generally can be said, that the maximum frequency deviation for the coordinated and the uncoordinated version show rather similar characteristics (cf. figure 12). In the case of standard MPC and passivity based MPC, the coordinated approach shows smaller frequency deviation than the uncoordinated approach, as expected and with inherent larger control input for coordinated version (cf. figure 16). Interestingly, with a horizon length of less than 38 steps, the Lyapunov based MPC setup exhibits a smaller maximum frequency deviation for the uncoordinated case. If a prediction horizon between 40 and 50 steps is applied, this swaps and the coordinated case shows a smaller maximum frequency deviation.
Naturally, as can be seen in figure 14, a smaller frequency deviation causes less power transmission between the areas.
Referring to the time for each optimization step (which excludes the optimizer setup time), the situation is similar to the one area setup, where CLF based and standard MPC needs comparably less time than passivity based MPC. This behavior seems to be fairly independent of the control topology. The coordinated CLF based MPC approach appears to be an exception, as the slope looks considerably steeper than for all other MPC approaches.
Vii Conclusion and Outlook
Viia Conclusion
In this paper a model predictive control based frequency control scheme for energy storage units was derived. The focus was on the incorporation of stability constraints, based on Lyapunov theory and the concept of passivity.
It has been shown that the concept of passivity as well as the concept of Control Lyapunov Functions (CLF) can easily be merged with the idea of MPC. The stability properties of the different MPC setups were derived, implemented and simulated. Furthermore the corresponding control performance was analyzed.
It was shown, that a standard MPC controller might evolve into unstable behavior. One of the reasons influencing this might be that energy storage charging level and the control input are penalized in the chosen setup. Augmenting this setup with stability constraints, such as in passivity based or CLF based MPC, stabilizes its behavior under otherwise equal conditions for all analyzed cases. To illustrate this in a study case, simulations were conducted for controlling a two area power system in a centralized and a decentralized setup. Using a longer prediction horizon might stabilize the conventional MPC controller, although this is usually not guaranteed.
Although the computational demands for passivity based MPC might be comparably higher than for standard MPC at equal prediction horizons, a considerably longer prediction horizon seems to be needed for achieving (not guaranteeing) stability in a standard MPC setup. This could lead to less computational costs for the passivity based MPC approach, e.g. a prediction horizon of only is needed in this setup with active passivity constraints to achieve acceptable control performance.
Regarding CLF based MPC, especially for the uncoordinated case, the time needed per optimization step is similar to the time needed within the standard MPC approach and results in guaranteed stability as in the passivity based approach.
Due to the fact, that asymptotic closed loop stability of the CLF based and passivity based MPC is guaranteed independently of the choice of and in the quadratic objective functional, a separation between performance and stability considerations is achieved to some degree. Specifically, the tuning of the corresponding MPC cost weights might be carried out with less effort to achieve acceptable control performance and stability.
ViiB Future Research
It might be promising to analyze how stability constraints (especially in a decentralized setup) would perform in a multinode power network. It is shown in [20], that for example in a threenode power network, the system may not be completely stable, even when the damping is made arbitrarily large – unlike in a twonode power network.
Furthermore, investigations on robustness properties of the proposed control schemes as well as including both active and reactive power transmission into the corresponding models might be of interest.
References
 [1] “Deutschland: Wind und Solar produzieren erstmals ï¿½ber 60 des Stroms.” www.eenews.ch/de/erneuerbare/article/26739/deutschlandwindundsolarproduzierenerstmalsueber6desstroms; Internationales Wirtschaftsforum Regenerative Energien (IWR), June 2013.
 [2] G. Andersson, Dynamics and Control of Electric Power Systems. EEHPower Systems Laboratory, ETH Zurich, Lecture Notes, 2012.
 [3] T. Rinke, A. Ulbig, S. Chatzivasileiadis, and G. Andersson, “Predictive control for realtime frequency control and inertia maximizing in power systems,” IEEE Transactions on Sustainable Energy, 2013.
 [4] P. Kundur, N. Balu, and M. Lauby, “Power system stability and control, EPRI power system engineering series, 1994.”
 [5] G. Andersson, Power System Analysis. EEHPower Systems Laboratory, ETH Zurich, Lecture Notes, 2012.
 [6] A. Ulbig, M. Galus, S. Chatzivasileiadis, and G. Andersson, “General frequency control with aggregated control reserve capacity from timevarying sources: The case of phevs,” in Bulk Power System Dynamics and Control (iREP)VIII (iREP), 2010 iREP Symposium, pp. 1–14, IEEE, 2010.
 [7] K. Heussen, S. Koch, A. Ulbig, and G. Andersson, “Unified systemlevel modeling of intermittent renewable energy sources and energy storage for power system operation,” Systems Journal, IEEE, vol. 6, no. 1, pp. 140–151, 2012.
 [8] T. Weißbach and E. Welfonder, “Improvement of the performance of scheduled stepwise power programme changes within the european power system,” in Proceedings of the 17th World Congress, The International Federation of Automatic Control (IFAC), Seoul, Korea, pp. 11972–11977, 2008.
 [9] M. Scherer, E. Iggland, A. Ritter, and G. Andersson, “Improved frequency bias factor sizing for noninteractive control,” CIGRE 2012, 2012.
 [10] M. Kurth and E. Welfonder, “Importance of the selfregulating effect within power systems,” IFAC Symposium on Power Plants and Power Systems Control, Kananaskis, Canada, pp. 345–352, 2006.
 [11] J. Carpentier, “’to be or not to be modern’ that is the question for automatic generation control (point of view of a utility engineer),” International Journal of Electrical Power & Energy Systems, vol. 7, no. 2, pp. 81–91, 1985.
 [12] D. Mayne, J. Rawlings, C. Rao, and P. Scokaert, “Constrained model predictive control: Stability and optimality,” Automatica, vol. 36, no. 6, pp. 789 – 814, 2000.
 [13] E. D. Sontag, “A universal construction of artstein’s theorem on nonlinear stabilization,” Systems & control letters, vol. 13, no. 2, pp. 117–123, 1989.
 [14] A. Jadbabaie, J. Yu, and J. Hauser, “Stabilizing receding horizon control of nonlinear systems: a control lyapunov function approach,” in American Control Conference, 1999. Proceedings of the 1999, vol. 3, pp. 1535–1539, IEEE, 1999.
 [15] J. A. Primbs, V. Nevistić, and J. C. Doyle, “Nonlinear optimal control: A control lyapunov function and receding horizon perspective,” Asian Journal of Control, vol. 1, no. 1, pp. 14–24, 1999.
 [16] T. Raff, C. Ebenbauer, and F. Allgöwer, “Passivitybased nonlinear model predictive control,” Assessment and Future Directions of Nonlinear Model Predictive Control, pp. 67–80, 2007.
 [17] R. Sepulchre, M. Janković, and P. Kokotović, Constructive nonlinear control. Communications and control engineering, Springer, 1997.
 [18] H. K. Khalil and J. Grizzle, Nonlinear systems, vol. 3. Prentice hall Upper Saddle River, 2002.
 [19] R. A. Freeman and P. V. Kokotovic, Robust nonlinear control design: statespace and Lyapunov techniques. Birkhauser, 2008.
 [20] A. Arapostathis and P. Varaiya, “Behaviour of threenode power networks,” International Journal of Electrical Power & Energy Systems, vol. 5, no. 1, pp. 22–30, 1983.