Patent ID: 12255462

DETAILED DESCRIPTION

When introducing elements of various embodiments of the present disclosure, the articles “a,” “an,” and “the” are intended to mean that there are one or more of the elements. The terms “comprising,” “including,” and “having” are intended to be inclusive and mean that there may be additional elements other than the listed elements. Additionally, it should be understood that references to “one embodiment” or “an embodiment” of the present disclosure are not intended to be interpreted as excluding the existence of additional embodiments that also incorporate the recited features. Furthermore, the phrase A “based on” B is intended to mean that A is at least partially based on B. Moreover, unless expressly stated otherwise, the term “or” is intended to be inclusive (e.g., logical OR) and not exclusive (e.g., logical XOR). In other words, the phrase “A or B” is intended to mean A, B, or both A and B.

In addition, several aspects of the embodiments described may be implemented as software modules or components. As used herein, a software module or component may include any type of computer instruction or computer-executable code located within a memory device and/or transmitted as electronic signals over a system bus or wired or wireless network. A software module or component may, for instance, include physical or logical blocks of computer instructions, which may be organized as a routine, program, object, component, data structure, or the like, and which performs a task or implements a particular data type.

In certain embodiments, a particular software module or component may include disparate instructions stored in different locations of a memory device, which together implement the described functionality of the module. Indeed, a module or component may include a single instruction or many instructions, and may be distributed over several different code segments, among different programs, and across several memory devices. Some embodiments may be practiced in a distributed computing environment where tasks are performed by a remote processing device linked through a communications network. In a distributed computing environment, software modules or components may be located in local and/or remote memory storage devices. In addition, data being tied or rendered together in a database record may be resident in the same memory device, or across several memory devices, and may be linked together in fields of a record in a database across a network.

Thus, embodiments may be provided as a computer program product including a tangible, non-transitory, computer-readable and/or machine-readable medium having stored thereon instructions that may be used to program a computer (or other electronic device) to perform processes described herein. For example, a non-transitory computer-readable medium may store instructions that, when executed by a processor of a computer system, cause the processor to perform certain methods disclosed herein. The non-transitory computer-readable medium may include, but is not limited to, hard drives, floppy diskettes, optical disks, compact disc read-only memories (CD-ROMs), digital versatile disc read-only memories (DVD-ROMs), read-only memories (ROMs), random access memories (RAMs), erasable programmable read-only memories (EPROMs), electrically erasable programmable read-only memories (EEPROMs), magnetic or optical cards, solid-state memory devices, or other types of machine-readable media suitable for storing electronic and/or processor executable instructions.

As discussed above, stability of a microgrid may be affected by numerous contingencies such as tripping of generation and load, faults in the power network or individual components, and the starting up of large induction machines. Conventional power delivery systems may use proportional-integral-derivative (PID) controllers and feedback loops to address stability issues occurring in a time period of a few seconds or less (e.g., less than three seconds, less than two seconds, and/or less than one second) (i.e., short-term stability). The PID controllers may measure parameters of the power delivery system and react to those measurements by adjusting excitation and turbine controls to change the output (e.g., voltage, frequency, and the like) of the power delivery system. However, the lower inertia and higher resistance-to-reactance ratios of conductors in a microgrid tend to worsen these instabilities.

Embodiments herein present control strategies to overcome the increased short-term instabilities of a microgrid. For example, embodiments presented herein may be based on a centralized model predictive control (MPC) for microgrid stability. The MPC predicts a change of a state of the power delivery system for a finite time period (e.g., a prediction period) based on a dynamic model. To do so, a centralized controller may obtain input data from various components of the power delivery system. The MPC may apply a trajectory linearization technique that separates nonlinear power system dynamics from adjustments of control input data. The MPC utilizes a mathematical model of dynamics of the power delivery system to predict development of future states and outputs of the system for the finite time period.

For example, the MPC may obtain system state data and perform calculations at a number Npof time steps which may occur every Tsseconds. In some embodiments, each time step k may be, for example, in a range of about 20 milliseconds and about 100 milliseconds, such as about 50 milliseconds or about 75 milliseconds. In some cases, the time step k may be less than 500 milliseconds. Based on the prediction, the MPC computes control inputs for the system over a control horizon having a number Ncof time steps. In some embodiments, an end of the control horizon may be the same as an end of the prediction horizon. In that case, a control input is computed for each step Ncof the control period. In some embodiments, the control period may be shorter than the prediction period. In that case, a number of variables, and thus a time to perform various calculations herein, may be reduced. In this way, a shorter sampling period (e.g., Ts) for each time step may reduce a number of variables and thus increase a speed of the calculations, thereby improving an efficiency and a performance of the MPC and the overall power delivery system. Alternatively, in some cases, the computational budget (e.g., computational resources to perform the calculations) may be maintained by reducing the sampling period with a reduced number of variables (e.g., time steps) without changing a length of the prediction horizon.

The power system model may be represented by the following nonlinear differential algebraic equations:
{dot over (x)}(t)=f(x(t),u(t),Iqd(t))  (1)
Iqd(t)=g(x(t))  (2)
y(t)=h(x(t),Iqd(t))  (3)
where x(t)∈nis a column vector of system states of the power delivery system, u(t)∈mis a vector of control inputs of the power delivery system, Iqd(t)∈2ris a column vector of qd-axis stator currents in individual rotor reference frames injected into a network of the power delivery system by r rotating machines, and y is a column vector of the outputs of the power delivery system which are tracked via the controller. The vector-valued function f represents the nonlinear differential equations associated with the system state dynamics. The function g represents algebraic network equations and reference frame transformations, and the function h represents algebraic system output equations. An assumption may be made that there is no direct feedthrough from input to output. These equations are the basis of the MPC and are discussed in more detail below. In some cases, an ordinary differential equation (ODE) solver may be used to obtain a continuous-time solution to the nonlinear differential-algebraic Equations (1)-(3).

The MPC may sample the power delivery system and initiates calculations every Tsseconds. At each time step k, the MPC may predict system states xk* and outputs yk* over the prediction horizon for a Nptime steps. A control input (e.g., an optimal control input, u(t)) may be computed over a control horizon of Nctime steps based on the predicted system states xk* and outputs yk*. It should be noted that an asterisk (e.g., xk*, yk*, and Iqd,k*) may be used to denote variables that are predicted by, for example, simulation.

Various delays may be caused by various operations executed by the power delivery system and/or the computations and communications of the MPC. These delays may be caused by, for example, measuring values such as voltage or current, performing calculations such as those discussed herein, updating network topology of the power delivery system, and/or communication between components of the power delivery system. In some cases, these delays may affect (e.g., increase and/or decrease) instabilities of the power delivery system. A worst case scenario of a sum of these delays may be denoted as τdelay.

The MPC may account for these delays using a delay period of Nd≥1 time steps, where NdTs≥τdelay. Thus, an optimization of the MPC that begins at time t0may include decision variables starting at t0+NdTsand may maintain preceding control inputs at an initial value (e.g., obtained from a previous optimization of the MPC at a previous time step, t0−Ts). To reduce a number of decision variable used by the MPC during computations, a length of the control horizon may be reduced (e.g., Nc<Np−Nd). The MPC may compute the control input u(t) for each time step Ncof the control horizon.

In some cases, the MPC may receive updates of the network topology from, for example, protective equipment in the network. For example, a protective relay may send a signal to the MPC to indicate that that relay cleared a fault by opening a line. The signal received by the MPC may indicate the open line. The MPC may account for a communication delay of the signal, as discussed above.

The MPC may begin calculating control inputs after the delay period. Thus, rather than the first calculation being performed at a first time step k, the MPC may perform the first calculation at time step k+Nd. In that case, the MPC may not calculate and/or change control inputs for time steps preceding the time step k+Nd. In this way, the MPC accounts for delays within the power delivery system and ensures control actions based on the control inputs are performed as needed to maintain or improve stability of the system and/or microgrid.

FIG.1is a schematic diagram of an electric power delivery system100, in accordance with an embodiment of the present disclosure. The electric power delivery system100may transmit power from generation to load using the depicted components. The electric power delivery system100may include may include generators110,112and loads114,116,118,120. Although the loads114,116,118,120are illustrated as motor loads, it should be understood that the system100may include and be coupled to various other types of loads, such as electrical equipment in an industrial plant. The generators110,112may generate power including active power and reactive power. Moreover, the loads114,116,118,120may consume the generated active power and the generated reactive power.

The generators110,112may be connected to a bus122via respective lines124,126(e.g., transmission lines and/or medium voltage lines) and via circuit breakers128,130. The circuit breakers128,130may be positioned on the respective lines124,126. Moreover, the loads114,116,118,120may be connected to the bus122via respective lines132,134,136,138and through circuit breakers140,142,144,146, and transformers148,150,152,154. The circuit breakers140,142,144,146and transformers148,150,152,154may be coupled to the respective lines132,134,136,138.

The electric power delivery system100may include a controller102coupled to one or more sensors108. In some embodiments, the controller102may include a model predictive control (MPC)104and may be a centralized controller for the system100. The one or more sensors108may obtain information related to the components of the electric power delivery system100, including the generators110,112, the loads114,116,118,120, the circuit breakers140,142,144,146, the transformers148,150,152,154, or any combination thereof. The controller102may monitor and control the power transmission of the electric power delivery system100.

As such, the controller102may receive one or more electrical measurements and/or measurements related to mechanical variables such as a generator rotor speed or angle, from the one or more sensors108. The sensors108may be proximate to the components of the power delivery system100and remote from the controller102. The controller102may also send control signals to the one or more components of the electric power delivery system100to control power transmission from generation to load (e.g., power flow). While two generators110,112and four loads114,116,118,120are depicted, it should be understood that the electric power delivery system may include more or fewer components than those shown.

In some embodiments, the controller102may open or close each of the circuit breakers128,130,140,142,144,146to control the power (active and/or reactive) generation and consumption of the electric power delivery system100. For example, during overcurrent conditions, the controller102may generate and transmit control signals to trip (open) one or more of the circuit breakers128,130,140,142,144,146to disconnect the generators110,112or the loads114,116,118,120. In specific embodiments, the controller102may generate and transmit such signals based on receiving inputs from the transformers148,150,152,154.

In some cases, the transformers148,150,152,154may be current transformers and may provide measurement of current flow through the respective lines132,134,136,138to the controller102. The controller102may use such current flow measurements to open or close one or more of the circuit breakers140,142,144,146associated with the loads114,116,118,120or circuit breakers128,130associated with the generators110,112. For example, the controller102may trip a circuit breaker (e.g., circuit breaker140,142,144,146) associated with a respective load (e.g., load114,116,118,120) with high current measurements to prevent overloading the respective load. Similarly, the controller102may also open or close one or more of the circuit breakers128,130associated with one or more generators110,112to prevent overloading different components of the electric power delivery system100.

Moreover, the controller102may receive measurements indicative of power (active and/or reactive) generated by the generators110,112and power delivered to the loads114,116,118,120. In some embodiments, the controller102may receive current measurements or voltage measurements of different components and translate such measurements to power measurements. Subsequently, the controller102may provide the control signals to respective controllers (e.g., governors, exciters) of one or more of the generators110,112or controllers of one or more loads114,116,118,120to adjust (e.g., increase and/or decrease) the power generation and the power consumption of the electric power delivery system100.

FIG.2is a flowchart depicting operations200for modeling an electric power delivery system and adjusting control inputs based on the model to improve or maintain stability of the electric power delivery system, according to an embodiment of the present disclosure. In some embodiments, the operations200may be performed by one or more processor and/or controllers, such as the controller102of the system100discussed with respect toFIG.1. It should be understood that, while the operations200are shown in a specific sequence, the operations200may be implemented in any suitable order, and at least some operations200may be skipped altogether.

The operations200begins at operation210where input data for the system100is obtained for a particular sampling period. The input data may be measured by one or more sensors108disposed proximate to an associated component of the system100. The input data may include a state of a system100, a voltage, a frequency, a rotor speed, a mechanical load torque, or any combination thereof. It should be noted that some of the input data may be computed based on other input data. For example, mechanical load torque may be calculated based on the rotor speed of an electric motor. The sampling period may occur at a regular interval and may have a duration that changes based on the input and/or output of the model predictive control (MPC). In some embodiments, the input data may include information associated with one or more loads of the associated microgrid.

At operation212, the controller102may obtain a nominal input trajectory and predict a state trajectory for (e.g., over) a prediction period (e.g., a prediction horizon). The nominal input trajectory over the prediction horizon may be a staircase signal represented as:
u*(t)=uk*  (4)
where t0+kTs≤t<t0+(k+1)Ts, ∀k∈{0, 1, . . . , Np−1}. The system state xk*, the output yk*, and qd currents Igd,k* may be obtained by evaluation of the continuous-time solution (e.g., obtained via an ODE solver) at time t0+kTsfor each discrete time step k of the prediction horizon.

At operation214, the MPC may obtain an input trajectory and predict a state trajectory for a prediction period (e.g., a prediction horizon). As discussed above, the continuous-time solution for each time step of the prediction period may be evaluated, via, for example, a simulation. That is, Jacobian matrices are evaluated at xk*, yk*, and Iqd,k*. Thus, Jacobian matrix elements may correspond to nonlinear terms recalculated for each time step k resulting in a linear state-space equation for the perturbed system dynamics for each time step k of the prediction horizon represented as:

ddt⁢δ⁢x⁡(t)=Ak⁢δ⁢x⁡(t)+Bk⁢δ⁢u⁡(t),(5)δ⁢y⁡(t)=Ck⁢δ⁢x⁡(t),(6)
where δ denotes a perturbation (e.g., a change) in the respective variable, and

Ak=∂f∂x+∂f∂Iqd⁢∂g∂x,Bk=∂f∂u,Ck=∂h∂x+∂h∂Iqd⁢∂g∂x.(7)

The MPC may discretize the continuous-time Equations (5)-(6) to obtain:
δxk=Φk-1δxk-1+Γk-1δuk-1,δx0=0.  (8)
δyk=Ckδxk,  (9)
where Φk=eAkTsand Γk=∫0TseAkηBkdη, for all k∈{1, . . . , Np}. By recursively applying Equations (8)-(9), perturbed system outputs may be written for each time step of the prediction period as a function of changes to the control input. Thus, the decision variables of the MPC are the perturbed control inputs:

δ⁢U=[δ⁢uNd⊤δ⁢uNd+1⊤…δ⁢uNd+Nc-1⊤]⊤∈ℝmNc.(10)

Manipulating the algebraic equations, a linear map between input and output perturbations may be obtained:

δ⁢Y=[δ⁢yNd+1⊤δ⁢yNd+2⊤…δ⁢yNp⊤]⊤=Z⁢δ⁢U,(11)
where δY∈p(Np−Nd), and Z is a block lower-triangular matrix. That is, the system dynamics may be linearized about each point (e.g., time step, k) of the state and output trajectories at operation214.

Linearization of the state trajectory and/or output trajectory enables the MPC to linearize the output trajectory at each point (e.g., time step) and reduce a complexity of the output (and state) trajectories. In this way, the MPC may improve an efficiency of the techniques described herein and improve a functioning of the MPC itself. Equation (8) represents how a change in either the control inputs (u) or system states (x) at the given time step k−1 impacts the change in system state (x) at the next time step (e.g., k). Equation (9) represents a change in system outputs (y) due to a change in the control inputs (u) or system states (x) for the given time step k.

At operation216, the MPC determines an input trajectory (u) for the prediction period based at least in part on the linearized output trajectory and state trajectory based on Equations (12)-(16) below. The input trajectory u may be determined by minimizing Equation (12) subject to the constraints of Equations (13)-(16). That is, after the output and state trajectories are linearized, the MPC may determine a change in the output trajectory (δy) in terms of a change in the input trajectory (du), for example, using Equation (11). The MPC may use this relationship in Equation (13) to solve Equations (12)-(16). Further, a first control input (u) that will affect the system output (y) occurs at Nd, which occurs sometime after the first time step (e.g., k=1).

The MPC may have two objectives: (1) maintain the output as close as possible to desired reference outputs and (2) return the system to a steady-state operating point as quickly as possible. These objections may be represented using the following objective function:

minδ⁢U,U,YJ=12⁢(Yref-Y)T⁢Q⁡(Yref-Y)+12⁢(H1⁢U)T⁢V1(H1⁢U)+12⁢(H2⁢δ⁢U)T⁢V2(H2⁢δ⁢U)(12)
such that
Y=Y*+ZδU,(13)
U=U*+δU(14)
ÃinU≤{tilde over (b)}in,  (15)
Umin≤U≤Umax,  (16)
where

H1=[-𝕀m𝕀m⋱⋱-𝕀m𝕀m],H2=[𝕀m0…0],(17)
and Y, Y*, Yref, U, U*, Umin, and Umaxare defined similarly to δY and δU of Equations (10) and (11) above.

The first term in Equation (12) (½(Yref−Y)TQ(Yref−Y)) may penalize deviations between system outputs and reference values, and thus is a key aspect in stabilizing the system during a transient. The second term (½(H1U)TV1(H1U)) and the third term (½(H2δU)TV2(H2δU)) in Equation (12) may be auxiliary terms and help to reduce oscillations in the decision variables once the system gets close to a steady-state. The second term of Equation (12) may penalize changes in control between consecutive time steps while the third term of Equation (12) is equal to

(1/2)⁢δ⁢uNdT⁢V2⁢δ⁢uNd.
A combined effect of the second and third terms at the steady-state is that uk=uNd* for all k≥Nd. The weight matrices Q, V1, and V2are diagonal and positive definite. In general, the elements of matrix Q should be chosen such that the first term dominates at the beginning of a transient. The constraints of Equations (15) and (16) may capture limits imposed by the control systems.

The objective function of Equation (12) may be re-written as a standard quadratic program equation:

minδ⁢UJ=12⁢δ⁢U⊤⁢P⁢δ⁢U+q⊤⁢δ⁢U+c(18)
such that
AinδU≤bin,  (19)
Umin−U*≤δU≤Umax−U*,(20)
where
P=ZTQZ+H1TV1H1+H2TV2H2,  (21)
g=−ZTQ(Yref−Y*)+H1TV1H1U*,(22)
c=½[(Yref−Y*)TQ(Yref−Y*)+(H1U*)TV1(H1U*)].  (23)

Equation (18) is convex, and thus has a unique solution, if and only if, P0. This solution can be verified by re-writing Equation (21) as
P=ZTQZ+HTVH,(24)
where HT=[H2TH1T] and V=blockdiag(V2,V1). Since V0 and H has a full column rank, HTVH0. Additionally, ZTQZ0 because Q0. Therefore, P0. Thus, the second and third terms of Equation (12) also ensure a unique solution to Equation (18), regardless of the rank of Z.

At operation218, the MPC may transmit control signals to components and/or equipment of the power delivery system based at least in part on the real-time control input values calculated at operation216. That is, the change in the input trajectory (e.g., (u) may be added to the nominal input trajectory u* to obtain a computed (e.g., “optimal”—acceptable, good, better, or best) input trajectory (e.g., u*+δu). The control signals transmitted to the power delivery system may be based at least in part on the computed input trajectory.

In some embodiments, the MPC may perform the operations200for each time step k (e.g., sampling period). That is, the MPC may obtain the input data, predict the nominal input trajectory, output trajectory and state trajectory, determine the updated input trajectory, determine real-time input data, and transmit the control signals based at least in part on the real-time input values, all within the time step. In some cases, the MPC may perform the operations200in a time period that is less than the time step k. In some embodiments, a remainder of the computed input trajectory for the current sampling period may be used as the initial predicted input trajectory for the next sampling period.

FIG.3is a graph250illustrating a trajectory linearization process of a model predictive control (MPC), according to an embodiment of the present disclosure. Specifically, input data u and output data y of a model predictive control (MPC) scheme and respective nominal input trajectory uk*262, predicted output trajectory yk*280, optimal input trajectory uk260, and optimal output trajectory yk282. A first portion252of the graph250illustrates the nominal input trajectory uk*262and corresponding optimal input trajectory uk260. The points266of the nominal input trajectory uk*262during the delay period256may represented by δuk=0, ∀k<Nd. The points268of the optimal input trajectory uk260after the control period258may be represented by uk=uNd+Nc−1, ∀k≥Nd+Nc.

A second portion254of the graph250illustrates the predicted output trajectory yk*280and the corresponding optimal output trajectory yk282. As shown, the graph of the250includes a delay period256ranging from time step k=0 to time step k=Nd286and a control period258(e.g., control horizon) ranging from the time step k=Nd286to the time step k=Nd+Nc288. A prediction period272(e.g., prediction horizon) for the delivery system100may have a duration from time step k=0 to the time step k=Np290.

As discussed above, the nominal input trajectory uk*262may be obtained by the MPC based on the input data and may be used to compute the optimal input trajectory uk260. The predicted output trajectory yk*280may be determined by solving the nonlinear ordinary differential Equations (1)-(3) above. After the output trajectory yk*280and the state trajectory x* are linearized as discussed with respect to operation214ofFIG.2, a change in the input trajectory δu264may be linearly mapped to a change in the output trajectory δy284as discussed above. In this way, the MPC may determine control inputs u to obtain the optimal output trajectory yk282.

As discussed above, two functions, among others, of the MPC may be to (1) maintain the output as close as possible to desired reference outputs and (2) return the system to a steady-state operating point as quickly as possible. The first function is illustrated in the second portion254ofFIG.3by the change in the output trajectory δy284showing that the predicted output trajectory yk*280is substantially similar to the corresponding optimal output trajectory yk282. Similarly, the second function of the MPC may be represented in the first portion252ofFIG.3by the points268showing that the optimal control input trajectory uk260returns to a steady state relatively quickly after the control period258.

FIG.4is a schematic diagram of an example microgrid300(e.g., electric power delivery system), according to an embodiment of the present disclosure. The example microgrid300may correspond to the electric power delivery system100ofFIG.1. As shown, the microgrid300includes a number of synchronous generators302,304,306,308coupled to busses310,312,314,316via components318, such as circuit breakers and lines320A,320B,322A,322B,324A, and324B. Each synchronous generator302,304,306,308consists of a synchronous machine with a field winding voltage controlled by a brushless exciter and mechanical torque controlled by either a gas or steam-driven turbine. As shown, a first synchronous generator302is a 33 MVA gas generator, a second synchronous generator304is a 60 MVA gas generator, a third synchronous generator308is a 153 MVA steam generator and a fourth synchronous generator306is a 19 MVA gas generator. It should be understood that the generators inFIG.4are merely examples and any type and/or combination of generators may be present in the microgrid300.

The busses310,312,314,316may be interconnected via one or more components318, such as transformers. As shown, a first bus310is coupled to the first synchronous generator302, a second bus312is coupled to the second synchronous generator304, a third bus316is coupled to the third synchronous generator308, and a fourth bus314is coupled to the fourth synchronous generator306. As shown, a load326,328,330,332is coupled to each bus310,312,314,316. In some embodiments, the loads326,328,330,332may be constant impedance loads and/or induction motors. As shown, the example microgrid300is a 13.8-kilovolt (kV) microgrid with a total generating capacity of 225 megawatts (MW).

To apply the MPC formulation above, the synchronous machines (e.g., generators) may be modeled using a “1.1 model” without stator transients. While a reduced order model is discussed below, it should be understood that a more detailed machine model(s) (e.g., a “2.2 model”) may be used in the MPC. The differential equations for the 1.1 model are:

Tdo′⁢adt⁢Eq′=-Eq′-(Xd-Xd′)⁢Id+Efd,(25)Tqo′⁢ddt⁢Ed′=-Ed′+(Xq-Xq′)⁢Iq,(26)ddt⁢δ=ωs(ω-1),(27)2⁢H⁢ddt⁢ω=Tm-Ed′⁢Id-Eq′⁢Iq-(Xq′-Xd′)⁢Id⁢Iq.(28)

All variables are per unit (pu), and ωsis the synchronous frequency. The rotor angle δ (not the perturbation operator) may be defined with respect to a common synchronously rotating reference frame. The variables Iqand Idare the stator currents injected into the network and are constrained by the algebraic network equations. The field winding voltage Efdand the mechanical torque Tmof Equations (25)-(28) may be inputs to the dynamic model of the synchronous machine. Linearizing the system dynamics fsm(xsm,usm,Iqd,sm) yields:

ddt⁢δ⁢xsm=∂fsm∂xsm⁢δ⁢xsm+∂fsm∂Iqd,sm⁢δ⁢Iqd,sm+∂fsm∂Tm⁢δ⁢Tm+∂fsm∂Efd⁢δ⁢Efd,(29)
where xsmT=[Eq′ Ed′ δ ω], usmT=[TmEfd], and Iqd,smT=[IqId]. The rotor field winding current Ifdis a linear function of states and stator currents. Thus,

δ⁢Ifd=∂Ifd∂xsm⁢δ⁢xsm+∂Ifd∂Iqd,sm⁢δ⁢Iqd,sm.(30)

FIG.5is a circuit diagram for an example voltage-behind-transient-reactance (VBR) circuit350, according to an embodiment of the present disclosure. The VBR circuit350may be representative of algebraic stator-side equations of the microgrid300, discussed with respect toFIG.4.FIG.6is a block diagram of an example brushless exciter380, according to an embodiment of the present disclosure. In some embodiments, the field winding voltage of the exciter alternator Efemay be controlled via a power amplifier circuit (e.g., with a digital excitation system), represented by a first-order transfer function.

The input voltage signal, Va, may be determined by the MPC. The limits on Vamay be accounted for in the controller using the constraints of Equation (20). The state Verepresents the voltage of the exciter 3-phase winding, which may be supplied to the rotating rectifier. The exciter dynamics are nonlinear due to saturation effects and may be dependent on the field winding current of the main alternator Ifd. The rectifier modes of operation may be captured by the following piecewise function:

fex(In)={1-33⁢In,0≤In≤3434-(In)2,34<In≤343⁢(1-In),34<In≤10,1<In.(31)

Linearizing the exciter dynamics fexc(xexc,uexc) yields

ddt⁢δ⁢xexc=∂fexc∂xexc⁢δ⁢xexc+∂fexc∂Va⁢δ⁢Va+∂fexc∂Ifd⁢δ⁢Ifd,(32)
where xexcT=[VeEfe] and uexcT=[VaIfd]. The exciter model outputs Efd=hexc(xexc,uexc)=VeFex. After algebraic manipulations,

δ⁢Efd=(Fex⋆-In⋆⁢∂Fex∂In)⁢δ⁢Ve+Ko⁢∂Fex∂In⁢δ⁢Ifd.(33)

Evaluating Equation (33) for the k-th interval of the prediction period272, it may be assumed that the exciter remains in the same mode of operation as predicted in the nonlinear simulation of Equations (1)-(3).

FIG.7is a block diagram of a gas and/or steam turbine model400used in a model predictive control (MPC), in accordance with an embodiment of the present disclosure. The model400may be representative of a simplified gas or steam turbine (with selection of proper parameters). The MPC may provide a reference fsr (e.g., fuel stroke demand) for a position of an actuator valve that controls a fuel flow Wf. The model400includes nonlinear slew rate limiters that limit a speed in which the actuator valve opens and/or closes. Conversion of the fuel into mechanical output power Pmby the turbine model400is modeled by a first-order transfer function. Mechanical torque Tmapplied to a shaft of the synchronous machine may be determined by dividing the mechanical power by the rotor speed. Linearizing the turbine dynamics fturb(xturb,uturb) yields

ddt⁢δ⁢xturb=∂fturb∂xturb⁢δ⁢xturb+∂fturb∂fsr⁢δ⁢fsr,(34)
where xturbT=[WfPm] and uturbT=[fsr ω]. The turbine model outputs Tm=hturb(xturb,uturb)=Pm/ω. Thus,

δ⁢Tm=Tm⋆(δ⁢PmPm⋆-δωω⋆).(35)

Limits on the fsr may be added to the MPC as constraints as in Equation (20). In some cases, the slew rate limiters in the turbine model may be ignored in computing the Jacobian matrices of Equation (34). However, impact of the slew rate limiters may be accounted for using the following constraints:
fsr−Wf≤TactropenNc(36)
−(fsr−Wf)≤−TactrcloseNc,  (37)
whereNcis a column vector of ones (1's). The vectors fsr and Wfmay contain the turbine input and fuel flows for each time step of the control period (e.g., fsr is a permutation of U). A procedure similar to that discussed with respect to Equations (4)-(11) to obtain
δWf=Lδfsr(38)
where L is a lower triangular matrix. Substitution yields trajectory-linearized inequality constraints of the form of Equation (15):

[𝕀Nc-LL-𝕀Nn]⁢δ⁢fsr≤[Tsct⁢ropen⁢𝕁Nc-(fsr⋆-Wf⋆)-Tact⁢rclose⁢𝕁Nc+(fsr⋆-Wf⋆)].(39)

FIG.8is a block diagram of an example synchronous generator model450used in the MPC, according to an embodiment of the present disclosure. The synchronous generator model450includes the exciter model380ofFIG.6, the turbine model400ofFIG.7, and a synchronous machine456. Linearized dynamics of the synchronous machine456of Equation (29), the exciter380of Equation (32), and the turbine400of Equation (34) may be combined in the general form of Equation (27), where:
xgen=[xsmT,xexcT,xturbT]T,ugen=[Va,fsr]T,  (40)
as shown inFIG.8.

For the constant impedance portion of the loads326,328,330,332ofFIG.4, an admittance of each load326,328,330,332may be added to a bus admittance matrix of the network. The induction machine models may consist of nonlinear differential algebraic equations similar to the synchronous machine models. An induction machine model may be used in the MPC formulation consisting four state variables (e.g., two for rotor flux linkages and two for mechanical dynamics) and an algebraic VBR circuit similar toFIG.5. The rotor windings of the induction machines may be assumed to be short-circuited and the mechanical load torque may be represented as a quadratic function of rotor speed. Similar to Equation (29) for the synchronous machine model, dynamics for the induction motors may be written in the form:
{dot over (x)}im=fim(xim,uim,Iqd,im)  (41)
where
xim=[Eq′,Ed′,δ,ω]T(42)
uim=Tm(43)
Iqd,im=[Iq,Id]T.  (44)

The power system network may be modeled without electromagnetic transients in the network components. Admittances of, for example, transmission lines, circuit breakers, transformers, etc., constant impedance loads, and VBR machine models in the microgrid300may be combined into a bus admittance matrix, {tilde over (Y)}. An impedance part of the VBR models350for the synchronous and induction machines may be added to the bus admittance matrix. For the synchronous generator VBR models350, the dependent voltage source may include sub-transient saliency terms that contain a q-axis stator current. For example, the qd-axis stator currents for the synchronous and induction machines may be expressed as a nonlinear algebraic function of the system states.

The stator currents for the r machines may be given by the nonlinear, complex algebraic equation:
(r−j{tilde over (Y)}S)T(x)Iq−jT(x)Id−{tilde over (Y)}T(x)E(x)=0,  (45)
where IqT=[Iq,1. . . Iq,r], IdT=[Id,1. . . Id,r], S∈r×ris a constant diagonal matrix that accounts for rotor saliency, T(x)=diag{[ejδ1. . . ejδr]} is the rotor reference frame transformation, and the VBR voltage sources are:
E(|x)=[Eq,1′−jEd,1′ . . . Eq,r′−jEd,r′]T.  (46)

Linearizing Equation (45) about the point (xk*,Iqd,k*) yields:
MqδIq+MdδId+Mxδx=0,  (47)
where Mq=(r−j{tilde over (Y)}S)T(xk*), Md=−jT(xk*), and

Mx=(𝕀r-j⁢Y~⁢S)⁢diag⁢{Iq,k⋆}⁢∂T∂x-j⁢diag⁢{Id,k⋆}⁢∂T∂x-Y~⁢T⁡(xk⋆)⁢∂E∂x-Y~⁢diag⁢{E⁡(xk⋆)}⁢∂T∂x.(48)
Equation 47 is split into real and imaginary parts, to obtain a sensitivity of the stator currents to system states, ∂g/∂x, which appears in the Jacobian matrices Akand Ckof Equation (7).

To ensure system stability, a proper selection of output values is needed. Presented herein are three types of short-term instability in AC microgrids: (1) transient (rotor angle) stability, (2) frequency stability, and (3) voltage stability. Transient instability may occur when rotors of synchronous machines lose synchronism. In one embodiment, to capture this phenomenon, differences in rotor speed of each generator may be examined with respect to a center-of-inertia (COI) speed (e.g., an inertia weighted average of rotor speed). If a rotor speed of a particular generator diverges from the COI speed, transient stability may have been lost.

Similarly, frequency instability may occur when a COI frequency of a generator diverges from a desired system frequency (e.g., 60 Hz) as defined by Equation (50) below. Voltage instability may occur when a magnitude of system bus voltage exceed a minimum limit or a maximum limit. These quantities may be captured in the MPC objective function for the microgrid300ofFIG.4using:

y=[ω^1…ω^4ωCOI❘"\[LeftBracketingBar]"V1❘"\[RightBracketingBar]"…❘"\[LeftBracketingBar]"V4❘"\[RightBracketingBar]"]T,(49)ω^g=ωg-ωCOI,ωCOI=∑g=14Hg⁢ωg∑g=14Hg,(50)❘"\[LeftBracketingBar]"Vi❘"\[RightBracketingBar]"=Vq,i2+Vd,i2=hvi(x,Iqd),(51)
where Hgis an inertia constant (e.g., in seconds) of generator g, Wo is a rotor speed of generator g, and Vq,iand Vd,iare qd-axis voltages at bus i. The VBR synchronous machine model may be used to express the voltage magnitude hvias a function of generator states and currents. The Jacobian matrix Ckin Equation (7) includes terms that are obtained by Equation (50) and linearizing Equation (51).

The reference output vector Yrefin Equation (12) may be set such that all terms associated with w of Equations (49)-(51) are 0 (e.g., over the prediction period), terms associated with ωCOImay be one (e.g., 1) per unit, and terms associated with |V| may be set to pre-fault voltage magnitudes determined by an optimal (e.g., acceptable, good, best) power flow (OPF) calculation. The weights in Q may be tuned to balance the desired importance of transient stability, frequency stability, voltage stability, and maintaining a steady-state operating point. The controlled inputs to the system are:
u=[fsr1Va,1. . . fsr4Va,4]T.  (52)

FIG.9is a block diagram of an example system model480for the microgrid300ofFIG.4, according to an embodiment of the present disclosure. The system model480may be generated by combining the dynamic models of the synchronous generators450ofFIG.8and induction machines (discussed with respect toFIG.8above) with the network model and system outputs. Despite the complexity of this system, Jacobian matrices such as those in Equation (7) may be computed for each MPC time step by making a few simplifications.

First, to address the piece-wise nature of the expression for Efdin the brushless exciter model380, it may be assumed that a mode of operation of the exciter does not change within the MPC time step for the linearization. Second, the slew rate limiters in the turbine model400were not included in computing the Jacobian matrices. However, the MPC may compensate for impacts caused by the slew rate limiters by using linear inequality constraints associated with the fuel stroke demand fsr to discourage the MPC from changing the fsr too quickly. Details of computing these Jacobian matrices are not included for brevity. However, it should be understood that any technique for computing Jacobian matrices, including the techniques above for computing the Jacobian matrices of Equation (7), could be used.

FIGS.10A-11Cillustrate example simulation results using an MPC as discussed herein for various faults in a microgrid, such as the microgrid300discussed with respect toFIG.3. As discussed below, the faults used for the simulation occur at various locations in the microgrid300.FIGS.10A-10Care based on data from all buses and generators of the microgrid300. As shown, all bus voltages are similar due to the relatively short line lengths of, for example, the lines320A,320B,322A,322B,324A, and324B. Voltage spikes may appear inFIGS.10A-10Cat the time (e.g., the instant, moment, etc.) when the fault is cleared due to snubbers in the circuit breaker models.

FIG.10Ais a graph500illustrating voltage magnitude profiles for all buses of the delivery system using the MPC discussed herein, according to an embodiment of the present disclosure.FIG.10Bis a graph512illustrating a center-of-inertia (COI) speed profile of the delivery system using the MPC discussed herein, according to an embodiment of the disclosure.FIG.10Cis a graph514illustrating a rotor speed deviation profile of the delivery system using the MPC, according to an embodiment of the present disclosure.

As shown, the graphs500,512, and514include a voltage magnitude profile for a maximum voltage magnitude502, a minimum voltage magnitude504, a median voltage magnitude506, a 95% prediction interval508, and a 75% prediction interval510. The MPC stabilizes system voltages and frequencies relatively quickly in each ofFIGS.10A-10C.

In some cases, operating conditions with heavily loaded induction motors near the fault may exhibit the worst performance, resulting in prolonged voltage sags and larger frequency deviations. The MPC performs well under all conditions, recovering system voltages within seconds and limiting frequency excursions to less than 1%.

FIGS.11A-11Care based on data from the third synchronous generator308. However, it should be understood that other generators in the microgrid may have similar responses. Further, as discussed above, the microgrid300is merely used as an example and it should be understood that different and/or additional components or a different microgrid topology may be used.

FIG.11Ais a graph520illustrating a center-of-inertia (COI) speed and a generator rotor speed in response to a fault on a corresponding bus of the delivery system using the MPC discussed herein, according to an embodiment of the present disclosure.FIG.11Bis a graph530illustrating a turbine input command and turbine fuel flow for the generator discussed with respect toFIG.11A, according to an embodiment of the present disclosure.FIG.11Cis a graph540illustrating an exciter input and bus voltage magnitude of the generator discussed with respect toFIGS.11A and111B, according to an embodiment of the present disclosure. For each ofFIGS.11A-11C, the graphs520,530,540illustrate responses for frequency (ω)522, center-of-inertia (COI) frequency (ωCOI)524, input voltage signal (Va)542, and the absolute value of the generator voltage V3544.

When the fault occurs (t=0.5 s), the voltage magnitude drops suddenly (FIG.11C) and the generator rotor speed initially falls, but then accelerates quickly (FIG.11A). At t=0.525 s, which is the first sampling point after occurrence of the fault, the MPC obtains knowledge of the fault from relays. The MPC updates the bus admittance matrix (e.g., the fault is modeled internally as a pure short-circuit to ground) and computes the optimal input. Due to delays, these control actions may be simulated as starting at the next sampling time, t=0.6 s, which is after the fault is cleared at t=0.5667 s. At t=0.6 s, one or more relays inform the MPC the fault is cleared and the admittance matrix is updated again to reflect the opening of faulted line (e.g., line324B). The MPC predicts the future drop in voltage magnitude and increase in rotor speed following the fault and responds by increasing the exciter input, Va, to its upper limit (FIG.11c) and decreasing fsr (FIG.11b). In general, the MPC decreases fsr from t=0.6 s to t=2.7 s when the COI frequency is high. However, during this interval, fsr is periodically increased as the rotor speed of the generator (e.g., the third generator308oscillates around the COI speed. The MPC respects turbine slew rate limits while adjusting the fsr.

The MPC presented herein improves response to a disturbance in a microgrid by utilizing power system models and real-time measurements of the microgrid from state observers and relays to provide a system-level perspective. Further, the MPC predicts a dynamic response of the system to disturbances rather than relying on a reactive feedback loop, tuning of which may be difficult and resource intensive. The MPC presented herein improves short-term frequency and voltage stability in microgrids. The trajectory linearization techniques discussed herein enables use of nonlinear component dynamic models in a real-time (or near real-time) implementation. In doing so, the trajectory linearization techniques provides an approach to solving complex, nonlinear power system dynamic models that is computationally efficient compared to other approaches.

As discussed above, the MPC control signals may be used to control inputs of a synchronous generator. In some embodiments, the MPC control signals may be used to control distributed energy resources (e.g., an inverter), inverter interfaced resources (e.g., power generated via solar, wind, batteries, or any combination thereof), control decisions such as load shedding, opening and/or closing circuit breakers and/or relays, and the like. To do so, these resources and/or components may be modeled via nonlinear differential algebraic equations similar to Equations (1)-(3) with state variables and inputs. Similar techniques regarding predictions, linearization, and determining real-time control input data may be used to achieve a desired output response.

It should be understood that the microgrids and associated components presented herein are merely examples and that the techniques presented herein may be applied to microgrids with more and/or fewer components and microgrids with additional and/or alternative components than those discussed herein. It should also be understood that the calculations and computations presented herein are merely example applications of the disclosed techniques and that fewer, additional, and/or alternative calculations and computations may be used to determine a change in an output of the components of the microgrid based on a change to the inputs of the components of the microgrid.

While specific embodiments and applications of the disclosure have been illustrated and described, it is to be understood that the disclosure is not limited to the precise configurations and components disclosed herein. For example, the systems and methods described herein may be applied to an industrial electric power delivery system or an electric power delivery system implemented in a boat or oil platform that may or may not include long-distance transmission of high-voltage power. Accordingly, many changes may be made to the details of the above-described embodiments without departing from the underlying principles of this disclosure. The scope of the present disclosure should, therefore, be determined only by the following claims.

Indeed, the embodiments set forth in the present disclosure may be susceptible to various modifications and alternative forms, specific embodiments have been shown by way of example in the drawings and have been described in detail herein. However, it may be understood that the disclosure is not intended to be limited to the particular forms disclosed. The disclosure is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the disclosure as defined by the following appended claims. In addition, the techniques presented and claimed herein are referenced and applied to material objects and concrete examples of a practical nature that demonstrably improve the present technical field and, as such, are not abstract, intangible or purely theoretical. Further, if any claims appended to the end of this specification contain one or more elements designated as “means for [perform]ing [a function] . . . ” or “step for [perform]ing [a function] . . . ”, it is intended that such elements are to be interpreted under 35 U.S.C. 112(f). For any claims containing elements designated in any other manner, however, it is intended that such elements are not to be interpreted under 35 U.S.C. 112(f).