# Stability analysis of a three-dimensional system of Topp model with diabetes

Z.S. Boxonov, U.A. Rozikov  U.Rozikova,b,c a V.I.Romanovskiy Institute of
Mathematics, 9, University str. 100174, Tashkent, Uzbekistan; b New Uzbekistan
University, 1, Movarounnahr str., 100000, Tashkent, Uzbekistan; c Faculty of
Mathematics, National University of Uzbekistan, 4, University str. 100174,
Tashkent, Uzbekistan<EMAIL_ADDRESS>Z.Boxonov, V.I.Romanovskiy Institute
of Mathematics, 9, University str. 100174, Tashkent, Uzbekistan.
<EMAIL_ADDRESS>
###### Abstract.

Mathematical models of glucose, insulin, and pancreatic $\beta$-cell mass
dynamics are essential for understanding the physiological basis of type 2
diabetes. This paper investigates the Topp model’s discrete-time dynamics to
represent these interactions. We perform a comprehensive analysis of the
system’s trajectory, examining both local and global behavior. First, we
establish the invariance of the positive trajectory and analyze the existence
of fixed points. Then, we conduct a complete stability analysis, determining
the local and global asymptotic stability of these fixed points. Finally,
numerical examples validate the effectiveness and applicability of our
theoretical findings. Additionally, we provide biological interpretations of
our results.

###### Key words and phrases:

glucose, insulin, beta-cell, fixed point, local stability, global behavior.

###### 2020 Mathematics Subject Classification:

39A12, 39A30, 92C50

##  1\. Introduction

To maintain normal human body functioning, it is essential to control the
level of glucose in the blood within the range of
$70-100mg/dl$111https://en.wikipedia.org/wiki/Blood-sugar-level. Insulin,
produced by beta-cells in the pancreas, facilitates the absorption of glucose
by cells and plays a critical role in blood glucose regulation. High blood
glucose levels trigger the release of insulin, which helps lower the
concentration to a healthy level. As blood glucose levels decrease, insulin
release gradually stops. This system, involving both insulin and glucagon,
helps maintain blood glucose balance and prevent diabetes-related
complications [2, 5, 10].

Mathematical models of glucose, insulin, and pancreatic beta-cell mass
dynamics are crucial for understanding the physiological basis of type 2
diabetes development. Traditionally, type 2 diabetes was thought to arise from
insulin insufficiency. The $GI\beta$ model developed by Topp and colleagues is
a prominent model for its progression. Topp’s model is foundational for
studying diabetes progression [2]. The Topp model is

$\left\\{\begin{array}[]{lll}\frac{dG}{dt}=g_{0}-g_{1}G-cGI,\\\\[8.53581pt]
\frac{dI}{dt}=\frac{s_{1}G^{2}}{s_{2}+G^{2}}\beta-kI,\\\\[8.53581pt]
\frac{d\beta}{dt}=(-d_{0}+r_{1}G-r_{2}G^{2})\beta,\end{array}\right.$ (1.1)

where $G$ $(mg/dl)$, $I$ $(\mu U/ml)$, $\beta$ $(mg)$ stand for the plasma
glucose concentration, insulin concentration and the mass of functional beta-
cells (preserving appropriate insulin production and secretion) at time $t$
(days), respectively. The parameter $g_{0}$ stands for the average rate of
glucose infusion per day (with meal ingestion as a main source), including the
hepatic glucose production. The term $g_{1}G$ represents insulin-independent
uptake of glucose, mainly by brain cells and nerve cells. In contrast, the
term $cGI$ depicts the insulin-dependent uptake of glucose, mostly by fat
cells and muscle cells in the human body. In particular, the coefficient $c$
$(ml/\mu U/day)$ stands for insulin sensitivity. The insulin secretion from
beta-cells is assumed to be triggered by increased glucose levels in the form
of the Hill function with coefficient 2, and the parameter $s_{1}$ represents
the secretory capacity per beta-cell. The insulin clearance rate is denoted by
$k$ (/day). The functional beta-cell mass is hypothesized to respond to
glucose with a pattern similar to a downward parabola: moderate amount of
glucose promotes the growth of beta-cells, while a high glucose level
exacerbates beta-cell apoptosis, resulting in the decrease of functional beta-
cell mass and $d_{0}$ $(d^{-1})$ is the death rate at zero glucose and $r_{1}$
$(mg^{-1}\ dl\ d^{-1})$, $r_{2}$ $(mg^{-2}\ dl^{2}\ d^{-1})$ are constants [2,
10].

In this paper (as in [6]-[9]) we study the discrete time dynamical systems
associated to the system (1.1). In the equations (1.1) we do the following
replacements:

$G=x,I=y,\beta=z.$

Define the operator $W:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ by

$\left\\{\begin{array}[]{lll}x^{\prime}=g_{0}-g_{1}x-cxy+x,\\\\[8.53581pt]
y^{\prime}=\frac{s_{1}x^{2}}{s_{2}+x^{2}}z-ky+y,\\\\[8.53581pt]
z^{\prime}=\left(-d_{0}+r_{1}x-r_{2}x^{2}\right)z+z,\end{array}\right.$ (1.2)

where all parameters in the model are positive. In the system (1.2),
$x^{\prime},y^{\prime},z^{\prime}$ means the next state relative to the
initial state $x,y,z$ respectively.

We call the partition into types hereditary if for each possible state
$\textbf{u}=(x,y,z)\in\mathbb{R}^{3}$ describing the current generation, the
state
$\textbf{u}^{\prime}=(x^{\prime},y^{\prime},z^{\prime})\in\mathbb{R}^{3}$ is
uniquely defined describing the next generation. This means that the
association $u\rightarrow u^{\prime}$ defines a map
$W:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ called the evolution operator [1],
[6].

The main problem for a given operator $W$ and arbitrarily initial point
$\textbf{u}^{(0)}=(x^{(0)},y^{(0)},z^{(0)})\in\mathbb{R}^{3}$ is to describe
the limit points of the trajectory $\\{\textbf{u}^{(m)}\\}_{m=0}^{\infty}$,
where

$\textbf{u}^{(m)}=W(\textbf{u}^{m-1})=W^{m}(\textbf{u}^{(0)}).$

##  2\. Local Stability Analysis of Fixed Points

Let $\mathbb{R}_{+}^{3}=\\{(x,y,z):x,y,z\in\mathbb{R},x\geq 0,y\geq 0,z\geq
0\\}.$

Note that the operator $W$ well defined on $\mathbb{R}^{3}$. But to define a
dynamical system of continuous operator as glucose, insulin, and pancreatic
beta-cell mass we assume $x\geq 0$, $y\geq 0$ and $z\geq 0$. Therefore, we
choose parameters of the operator $W$ to guarantee that it maps
$\mathbb{R}^{3}_{+}$ to itself.

We denote

$\Omega=\left\\{(x,y,z)\in\mathbb{R}_{+}^{3}:g_{0}\leq x\leq A,0\leq y\leq
B,0\leq z\leq C\right\\},$

where

$A=\frac{r_{1}+\sqrt{r_{1}^{2}+4r_{2}(1-d_{0})}}{2r_{2}},\
B=\frac{1-g_{1}}{c},\ C=\frac{(1-g_{1})k}{s_{1}c}.$

######  Lemma 1.

If

$0<g_{1}<1,\ 0<k\leq 1,\ 0<d_{0}\leq 1,\ g_{0}/g_{1}\leq A,\ r_{1}^{2}\leq
4r_{2}d_{0}$ (2.1)

then the operator (1.2) maps the set $\Omega$ to itself.

###### Proof.

Let for any $(x,y,z)\in\Omega$, i.e., $g_{0}\leq x\leq A,0\leq y\leq B,0\leq
z\leq C.$ Then

$x^{\prime}=g_{0}+x(1-g_{1}-cy)\geq g_{0}+x(1-g_{1}-c\cdot
B)=g_{0}+x(1-g_{1}-c\cdot\frac{1-g_{1}}{c})=g_{0},$
$x^{\prime}=g_{0}+x(1-g_{1}-cy)\leq g_{0}+x(1-g_{1})\leq
g_{0}+A(1-g_{1})=A-g_{1}(A-g_{0}/g_{1})\leq A.$

Thus $g_{0}\leq x^{\prime}\leq A.$

Since $0<k\leq 1$, it is clear that $y^{\prime}\geq 0.$ We will now show that
$y^{\prime}\leq B.$

$y^{\prime}=\frac{s_{1}x^{2}}{s_{2}+x^{2}}z-ky+y\leq s_{1}\cdot
C-ky+B=s_{1}\cdot\frac{(1-g_{1})k}{s_{1}c}-ky+B=B-k(B-y)\leq B.$

It follows that $0\leq y^{\prime}\leq B.$

It can be seen that the number $A$ is a positive solution to the equation
$1-d_{0}+r_{1}x-r_{2}x^{2}=0$. When $x$ ranges from $0$ to $A$, the expression
$1-d_{0}+r_{1}x-r_{2}x^{2}$ is always non-negative. Since
$[g_{0},A]\subset[0,A]$, then $z^{\prime}\geq 0.$

$z^{\prime}=\left(1-d_{0}+r_{1}x-r_{2}x^{2}\right)z\leq\left(1-d_{0}+\frac{r_{1}^{2}}{2r_{2}}-\frac{r_{1}^{2}}{4r_{2}}\right)C=\left(1-\frac{4r_{2}d_{0}-r_{1}^{2}}{4r_{2}}\right)C\leq
C.$

It follows that $z^{\prime}$ also changes from 0 to $C$. Thus
$(x^{\prime},y^{\prime},z^{\prime})\in\Omega.$ ∎

###  2.1. Fixed points.

First, we discuss the existence of the fixed points. A point $u\in\Omega$ is
called a fixed point of $W$ if $W(u)=u$.

######  Proposition 1.

The fixed points for (1.2) are as follows:

  * •

If $r_{1}^{2}<4r_{2}d_{0}$ then the operator (1.2) has a unique fixed point
$u^{*}=\left(\frac{g_{0}}{g_{1}},0,0\right).$

  * •

If $r_{1}^{2}=4r_{2}d_{0}$,
$\frac{g_{1}r_{1}}{2r_{2}}<g_{0}<\frac{r_{1}}{2r_{2}}$,
$s_{2}r_{2}(2r_{2}g_{0}-g_{1}r_{1})\leq d_{0}(r_{1}-2r_{2}g_{0})$ then mapping
(1.2) has two fixed points with

$u^{*}_{1}=\left(\frac{g_{0}}{g_{1}},0,0\right),\ \
u^{*}_{2}=\left(x^{*},y^{*},z^{*}\right),$

where

$x^{*}=\frac{r_{1}}{2r_{2}},\ y^{*}=\frac{g_{0}-g_{1}x^{*}}{cx^{*}},\
z^{*}=\frac{ky^{*}(s_{2}+x^{*}{{}^{2}})}{s_{1}x^{*}{{}^{2}}}.$ (2.2)

###### Proof.

The equation $W(u)=u$ is the following system

$\left\\{\begin{array}[]{lll}x=g_{0}+x(1-g_{1}-cy),\\\\[8.53581pt]
y=\frac{s_{1}x^{2}}{s_{2}+x^{2}}z+(1-k)y,\\\\[8.53581pt]
z=\left(1-d_{0}+r_{1}x-r_{2}x^{2}\right)z.\end{array}\right.$ (2.3)

The third equation of system (2.3) shows that either $z=0$ or
$r_{2}x^{2}-r_{1}x+d_{0}=0.$ It is easy to see that
$x=\frac{g_{0}}{g_{1}},y=0,z=0$ and
$x^{*}=\frac{r_{1}}{2r_{2}},y^{*}=\frac{g_{0}-g_{1}x^{*}}{cx^{*}},z^{*}=\frac{ky^{*}(s_{2}+x^{*}{{}^{2}})}{s_{1}x^{*}{{}^{2}}}$
(under conditions (2.1)) are solution to (2.3). Now let’s find additional
conditions for the parameters so that the point $(x^{*},y^{*},z^{*})$ belongs
to the set $\Omega$.

If $r_{1}^{2}<4r_{2}d_{0}$ then $x^{*}\notin[g_{0},A],$ if
$r_{1}^{2}=4r_{2}d_{0}$ and $2r_{2}g_{0}<r_{1}$ then $x^{*}\in[g_{0},A].$ It
follows that $y\leq B.$ Indeed,

$y^{*}=\frac{g_{0}-g_{1}x^{*}}{cx^{*}}\leq\frac{x^{*}-g_{1}x^{*}}{cx^{*}}=\frac{1-g_{1}}{c}=B.$

If $g_{1}r_{1}<2r_{2}g_{0}$ then $y^{*}\geq 0.$, i.e., $y^{*}\in[0,B].$ Since
$x$ and $y$ are non-negative, it follows that $z$ is also non-negative. We
solve the inequality $z^{*}\leq C$ and form the condition
$s_{2}r_{2}(2r_{2}g_{0}-g_{1}r_{1})\leq d_{0}(r_{1}-2r_{2}g_{0})$ for the
parameters.

So, if the parameters satisfy the following conditions

$r_{1}^{2}=4r_{2}d_{0},\
\frac{g_{1}r_{1}}{2r_{2}}<g_{0}<\frac{r_{1}}{2r_{2}},\
s_{2}\leq\frac{d_{0}(r_{1}-2r_{2}g_{0})}{r_{2}(2r_{2}g_{0}-g_{1}r_{1})},$
(2.4)

then the point $(x^{*},y^{*},z^{*})$ belongs to the set $\Omega$. ∎

###  2.2. Types of the fixed points

Now we shall examine the type of the fixed points.

######  Definition 1.

(see [3]) A fixed point $u^{*}$ of the operator $W$ is called hyperbolic if
its Jacobian $J$ at $u^{*}$ has no eigenvalues on the unit circle.

######  Definition 2.

(see [3]) A hyperbolic fixed point $u^{*}$ called:

1) attracting if all the eigenvalues of the Jacobi matrix $J(u^{*})$ are less
than 1 in absolute value;

2) repelling if all the eigenvalues of the Jacobi matrix $J(u^{*})$ are
greater than 1 in absolute value;

3) a saddle otherwise.

Before analyzing the fixed points we give the following useful lemma ([4]).

######  Lemma 2.

Let $F(\lambda)=\lambda^{2}+B^{*}\lambda+C^{*},$ where $B^{*}$ and $C^{*}$ are
two real constants. Suppose $\lambda_{1}$ and $\lambda_{2}$ are two roots of
$F(\lambda)=0$. Then the following statements hold.

  1. (i)

If $F(1)>0$ then

  2. (i.1)

$|\lambda_{1}|<1$ and $|\lambda_{2}|<1$ if and only if $F(-1)>0$ and
$C^{*}<1;$

  3. (i.2)

$\lambda_{1}=-1$ and $\lambda_{2}\neq-1$ if and only if $F(-1)=0$ and
$B^{*}\neq 2;$

  4. (i.3)

$|\lambda_{1}|<1$ and $|\lambda_{2}|>1$ if and only if $F(-1)<0;$

  5. (i.4)

$|\lambda_{1}|>1$ and $|\lambda_{2}|>1$ if and only if $F(-1)>0$ and
$C^{*}>1;$

  6. (i.5)

$\lambda_{1}$ and $\lambda_{2}$ are a pair of conjugate complex roots and
$|\lambda_{1}|=|\lambda_{2}|=1$ if only if $-2<B^{*}<2$ and $C^{*}=1;$

  7. (i.6)

$\lambda_{1}=\lambda_{2}=-1$ if only if $F(-1)=0$ and $B=2.$

  8. (ii)

If $F(1)=0,$ namely, 1 is one root of $F(\lambda)=0,$ then the other root
$\lambda$ satisfies $|\lambda|=(<,>)1$ if and only if $|C^{*}|=(<,>)1.$

  9. (iii)

If $F(1)<0$ then $F(\lambda)=0$ has one root lying in $(1;\infty).$ Moreover,

  10. (iii.1)

the other root $\lambda$ satisfies $\lambda<(=)-1$ if and only if
$F(-1)<(=)0;$

  11. (iii.2)

the other root $\lambda$ satisfies $-1<\lambda<1$ if and only if $F(-1)>0.$

To find the type of a fixed point of the operator (1.2) we write the Jacobi
matrix:

$J(u)=J_{W}=\left(\begin{array}[]{ccc}1-g_{1}-cy&-cx&0\\\
\frac{2s_{1}s_{2}x}{(s_{2}+x^{2})^{2}}z&1-k&\frac{s_{1}x^{2}}{s_{2}+x^{2}}\\\
(r_{1}-2r_{2}x)z&0&1-d_{0}+r_{1}x-r_{2}x^{2}\\\ \end{array}\right).$

The eigenvalues of the Jacobi matrix at the fixed point
$u_{1}^{*}=\left(\frac{g_{0}}{g_{1}},0,0\right)$ are as follows

$\lambda_{1}=1-g_{1},\ \lambda_{2}=1-k,\
\lambda_{3}=1-d_{0}+\frac{r_{1}g_{0}}{g_{1}}-\frac{r_{2}g_{0}^{2}}{g^{2}_{1}}.$

By (2.1) we have $0<\lambda_{1}<1,0\leq\lambda_{2}<1$. We write $\lambda_{3}$
as follows.

$\lambda_{3}=1-d_{0}+\frac{r_{1}g_{0}}{g_{1}}-\frac{r_{2}g_{0}^{2}}{g^{2}_{1}}=1-r_{2}\left(\frac{g_{0}}{g_{1}}-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r^{2}_{1}-4r_{2}d_{0}}{4r_{2}}.$

If $r^{2}_{1}\leq 4r_{2}d_{0}$ then $\lambda_{3}<1$.

We calculate eigenvalues of Jacobian matrix at the fixed point
$u_{2}^{*}=\left(x^{*},y^{*},z^{*}\right)$. The characteristic equation is

$(\lambda-1)F(\lambda)=(\lambda^{2}+B^{*}\lambda+C^{*})(\lambda-1),$ (2.5)

where

$B^{*}=k+\frac{g_{0}}{x^{*}}-2,$

$C^{*}=\left(1-\frac{g_{0}}{x^{*}}\right)(1-k)+\frac{2s_{2}kcy^{*}}{s_{2}+x^{*}{{}^{2}}}.$

From (2.5) it is clear that one of the eigenvalues is equal to one. Let us
determine the sign of $F(1)$, $F(-1)$ and $C^{*}-1$ under conditions (2.1),
(2.4).

$F(1)=1+B^{*}+C^{*}=1+k+\frac{g_{0}}{x^{*}}-2+\left(1-\frac{g_{0}}{x^{*}}\right)(1-k)+\frac{2s_{2}kcy^{*}}{s_{2}+x^{*}{{}^{2}}}=\frac{g_{0}k}{x^{*}}+\frac{2s_{2}kcy^{*}}{s_{2}+x^{*}{{}^{2}}}>0,$

$F(-1)=1-B^{*}+C^{*}=1-k-\frac{g_{0}}{x^{*}}+2+\left(1-\frac{g_{0}}{x^{*}}\right)(1-k)+\frac{2s_{2}kcy^{*}}{s_{2}+x^{*}{{}^{2}}}=\left(2-\frac{g_{0}}{x^{*}}\right)(2-k)+\frac{2s_{2}kcy^{*}}{s_{2}+x^{*}{{}^{2}}}>0,$

$C^{*}-1=\left(1-\frac{g_{0}}{x^{*}}\right)(1-k)+\frac{2s_{2}kcy^{*}}{s_{2}+x^{*}{{}^{2}}}-1=\frac{2s_{2}k(g_{0}-g_{1}x^{*})}{x^{*}(s_{2}+x^{*}{{}^{2}})}-\frac{g_{0}}{x^{*}}\left(1-k\right)-k=$

$=\frac{1}{x^{*}(s_{2}+x^{*}{{}^{2}})}\cdot\left(2s_{2}k(g_{0}-g_{1}x^{*})-g_{0}(1-k)(s_{2}+x^{*}{{}^{2}})-kx^{*}(s_{2}+x^{*}{{}^{2}})\right)\leq$

$\leq\frac{k}{x^{*}(s_{2}+x^{*}{{}^{2}})}\cdot\left(2s_{2}(g_{0}-g_{1}x^{*})-x^{*}(s_{2}+x^{*}{{}^{2}})\right)=$

$=\frac{k}{x^{*}(s_{2}+x^{*}{{}^{2}})}\cdot\left(s_{2}(2(g_{0}-g_{1}x^{*})-x^{*})-x^{*}{{}^{3}}\right)\leq$

$\leq\frac{k}{x^{*}(s_{2}+x^{*}{{}^{2}})}\cdot\left(\frac{x^{*}{{}^{2}}(x^{*}-g_{0})}{g_{0}-g_{1}x^{*}}\cdot\left(2(g_{0}-g_{1}x^{*})-x^{*}\right)-x^{*}{{}^{2}}\right)=$

$=\frac{kx^{*}}{(s_{2}+x^{*}{{}^{2}})(g_{0}-g_{1}x^{*})}\cdot\left(-(1+g_{1})x^{*}{{}^{2}}+2g_{0}(1+g_{1})x^{*}-2g_{0}^{2}\right)=$

$=-\frac{k(1+g_{1})x^{*}}{(s_{2}+x^{*}{{}^{2}})(g_{0}-g_{1}x^{*})}\cdot\left((x^{*}-g_{0})^{2}+\frac{(1-g_{1})g_{0}^{2}}{1+g_{1}}\right)<0.$

According to item (i.1) of Lemma 2, $|\lambda_{2,3}|<1.$

Thus the type of fixed points the following theorem holds.

######  Theorem 1.

The type of the fixed points for (1.2) are as follows:

  * i)

if $r_{1}^{2}<4r_{2}d_{0}$, then the unique fixed point $u^{*}$ of the
operator (1.2) is the attracting.

  * ii)

if $r_{1}^{2}=4r_{2}d_{0}$,
$\frac{g_{1}r_{1}}{2r_{2}}<g_{0}<\frac{r_{1}}{2r_{2}}$,
$s_{2}r_{2}(2r_{2}g_{0}-g_{1}r_{1})\leq d_{0}(r_{1}-2r_{2}g_{0})$ then the
operator (1.2) has two fixed points $u_{1}^{*},u_{2}^{*}$ and the point
$u_{1}^{*}$ is attracting, the point $u_{2}^{*}$ is non-hyperbolic (but,
$u_{2}^{*}$ is semi-attracting 222means that two eigenvalues are less than 1
in absolute value.).

##  3\. Periodic points.

A point u in $\Omega$ is called periodic point of $W$ if there exists
$p\in\mathbb{N}$ so that $W^{p}(\textbf{u})=\textbf{u}$. The smallest positive
integer $p$ satisfying $W^{p}(\textbf{u})=\textbf{u}$ is called the prime
period or least period of the point $\textbf{u}.$

######  Theorem 2.

For $p\geq 2$ the operator (1.2) does not have any $p$-periodic point in the
set $\Omega.$

###### Proof.

Let’s consider the following system.

$\left\\{\begin{array}[]{lll}x=x^{(p)}=g_{0}+x^{(p-1)}\left(1-g_{1}-cy^{(p-1)}\right),\\\\[8.53581pt]
y=y^{(p)}=\frac{s_{1}\left(x^{(p-1)}\right){{}^{2}}}{s_{2}+\left(x^{(p-1)}\right){{}^{2}}}z^{(p-1)}+(1-k)y^{(p-1)},\\\\[8.53581pt]
z=z^{(p)}=\left(1-d_{0}+r_{1}x^{(p-1)}-r_{2}\left(x^{(p-1)}\right){{}^{2}}\right)z^{(p-1)}.\end{array}\right.$
(3.1)

From the third equation of system (3.1) we have

$z=z^{(p-1)}\left(1-r_{2}\left(x^{(p-1)}-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)=z^{(p-2)}\left(1-r_{2}\left(x^{(p-1)}-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)\left(1-r_{2}\left(x^{(p-2)}-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)=...=z\left(1-r_{2}\left(x^{(p-1)}-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)\left(1-r_{2}\left(x^{(p-2)}-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)\cdot...\cdot\left(1-r_{2}\left(x-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right).$

Since $z\neq 0$, we have

$\begin{array}[]{ll}\left(1-r_{2}\left(x^{(p-1)}-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)\left(1-r_{2}\left(x^{(p-2)}-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)\cdot...\\\
...\cdot\left(1-r_{2}\left(x-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)=1\end{array}$
(3.2)

Under the condition $r_{1}^{2}<4r_{2}d_{0}$, the left side of equation (3.2)
will always be less than 1. In this case, the equation will not have roots,
which means that operator (1.2) does not have a $p$-periodic point in the set
$\Omega$.

Let $r_{1}^{2}=4r_{2}d_{0}.$ Then, for equation (3.2) to have a root, it is
necessary and sufficient to satisfy the relations

$x^{(p-1)}=x^{(p-2)}=...=x=\frac{r_{1}}{2r_{2}}.$

These relationships hold only at a fixed point $\frac{r_{1}}{2r_{2}}.$ So, in
this case, operator (1.2) does not have $p$-periodic points in the set
$\Omega$. ∎

##  4\. Global Behavior

In this section for any initial point
$\left(x^{(0)},y^{(0)},z^{(0)}\right)\in\Omega$ we investigate behavior of the
trajectories
$\left(x^{(n)},y^{(n)},z^{(n)}\right)=W^{n}\left(x^{(0)},y^{(0)},z^{(0)}\right),n\geq
1.$

We have

$\begin{array}[]{lll}x^{(n)}=g_{0}+x^{(n-1)}\left(1-g_{1}-cy^{(n-1)}\right),\\\\[8.53581pt]
y^{(n)}=\frac{s_{1}\left(x^{(n-1)}\right){{}^{2}}}{s_{2}+\left(x^{(n-1)}\right){{}^{2}}}z^{(n-1)}+(1-k)y^{(n-1)},\\\\[8.53581pt]
z^{(n)}=\left(1-d_{0}+r_{1}x^{(n-1)}-r_{2}\left(x^{(n-1)}\right){{}^{2}}\right)z^{(n-1)}.\end{array}$
(4.1)

######  Lemma 3.

Assume that (2.1) holds. For sequences $x^{(n)}$, $y^{(n)}$ and $z^{(n)}$, the
following property holds:

  * (i)

if $r_{1}^{2}<4r_{2}d_{0}$, then the sequence $z^{(n)}$ is monotonically
decreasing and converges to zero;

  * (ii)

if the sequence $z^{(n)}$ converges to zero, then the sequence $y^{(n)}$ also
has a limit and converges to zero;

  * (iii)

if the sequence $y^{(n)}$ converges to zero, then the sequence $x^{(n)}$ also
has a limit and converges to $\frac{g_{0}}{g_{1}}$.

###### Proof.

First, we prove the assertion (i). From third equation of system (4.1) we get

$\frac{z^{(n)}}{z^{(n-1)}}=1-d_{0}+r_{1}x^{(n-1)}-r_{2}\left(x^{(n-1)}\right){{}^{2}}=1-r_{2}\left(x^{(n-1)}-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}<1$.

This implies that $z^{(n)}$ is a decreasing sequence.

$z^{(n)}=\left(1-d_{0}+r_{1}x^{(n-1)}-r_{2}\left(x^{(n-1)}\right){{}^{2}}\right)z^{(n-1)}=\left(1-r_{2}\left(x^{(n-1)}-\frac{r_{1}}{2r_{2}}\right)^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)z^{(n-1)}<\left(1+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)z^{(n-1)}<\left(1+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)^{2}z^{(n-2)}<...<\left(1+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)^{n}z^{(0)}.$

Thus $0\leq
z^{(n)}<\left(1+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)^{n}z^{(0)}.$
Consequently

$\lim\limits_{n\to\infty}z^{(n)}=0.$

Let’s prove the assertion (ii). By analyzing the second equation of (4.1) and
considering the boundedness of the sequence $x^{(n)}$, we have

$\frac{s_{1}g_{0}^{2}}{s_{2}+g_{0}^{2}}z^{(n-1)}\leq
y^{(n)}-(1-k)y^{(n-1)}\leq\frac{s_{1}A^{2}}{s_{2}+A^{2}}z^{(n-1)}.$ (4.2)

Since the sequence $z^{(n)}$ converges to zero, it follows that the sequence
$y^{(n)}-(1-k)y^{(n-1)}$ also converges to zero.

Since $y^{(n)}$ is bounded, it has a well-defined upper limit, denoted by
$\alpha$. Then there must exist a subsequence $y^{(n_{j})}$ that converges to
$\alpha$. Furthermore, the fact that $y^{(n)}-(1-k)y^{(n-1)}$ converges to
zero implies another subsequence $y^{(n_{j}-1)}$ approaches
$\frac{\alpha}{1-k}$. Since $k$ is between 0 and 1, $\frac{\alpha}{1-k}$ must
be greater than or equal to $\alpha$. However, since $\alpha$ is an upper
limit, it follows that $\frac{\alpha}{1-k}\leq\alpha.$ Therefore,
$\frac{\alpha}{1-k}$ can only equal $\alpha$, which implies $\alpha=0$. We
conclude that the sequence $y^{(n)}$ has a limit, which is equal to zero.

(iii). Using the first equation of (4.1) and the fact that $x^{(n)}$ is
bounded, we can conclude that

$g_{0}-cAy^{(n-1)}\leq x^{(n)}-(1-g_{1})x^{(n-1)}\leq g_{0}-cg_{0}y^{(n-1)}.$
(4.3)

Since $y^{(n)}$ converges to zero, it follows from (4.3) that

$\lim\limits_{n\to\infty}\left(x^{(n)}-(1-g_{1})x^{(n-1)}\right)=g_{0}.$

Since $x^{(n)}$ is bounded, it possesses upper and lower limits, denoted by
$\alpha$ and $\beta$, respectively. To arrive at a contradiction, suppose
$x^{(n)}$ has no limit (i.e.,$\alpha\neq\beta$). This would necessitate the
existence of two subsequences: $x^{(n_{i})}$ converging to $\alpha$ and
$x^{(n_{j})}$ converging to $\beta$.

$\lim\limits_{i\to\infty}\left(x^{(n_{i})}-(1-g_{1})x^{(n_{i}-1)}\right)=g_{0},\
\Rightarrow\ \lim\limits_{i\to\infty}x^{(n_{i}-1)}=\frac{\alpha-
g_{0}}{1-g_{1}}\leq\alpha,\ \Rightarrow\ \alpha\leq\frac{g_{0}}{g_{1}},$
$\lim\limits_{i\to\infty}\left(x^{(n_{i}+1)}-(1-g_{1})x^{(n_{i})}\right)=g_{0},\
\Rightarrow\
\lim\limits_{i\to\infty}x^{(n_{i}+1)}=g_{0}+(1-g_{1})\alpha\leq\alpha,\
\Rightarrow\ \alpha\geq\frac{g_{0}}{g_{1}},$
$\lim\limits_{j\to\infty}\left(x^{(n_{j})}-(1-g_{1})x^{(n_{j}-1)}\right)=g_{0},\
\Rightarrow\ \lim\limits_{j\to\infty}x^{(n_{j}-1)}=\frac{\beta-
g_{0}}{1-g_{1}}\geq\beta,\ \Rightarrow\ \beta\geq\frac{g_{0}}{g_{1}},$
$\lim\limits_{j\to\infty}\left(x^{(n_{j}+1)}-(1-g_{1})x^{(n_{j})}\right)=g_{0},\
\Rightarrow\
\lim\limits_{j\to\infty}x^{(n_{j}+1)}=g_{0}+(1-g_{1})\beta\geq\beta,\
\Rightarrow\ \beta\leq\frac{g_{0}}{g_{1}},$

Therefore, we conclude that $\alpha=\beta=\frac{g_{0}}{g_{1}}$ i.e.,
$\lim\limits_{n\to\infty}x^{(n)}=\frac{g_{0}}{g_{1}}$. ∎

######  Theorem 3.

Assume that (2.1) holds and let $r_{1}^{2}<4r_{2}d_{0}$. Then the trajectory
converges to fixed point $u^{*}$ for any initial point
$\left(x^{(0)},y^{(0)},z^{(0)}\right)\in\Omega$, i.e.,

$\lim\limits_{n\to\infty}W^{n}\left(x^{(0)},y^{(0)},z^{(0)}\right)=\left(\frac{g_{0}}{g_{1}},0,0\right).$

###### Proof.

The proof follows from Lemma 3. ∎

A set $S$ is called invariant with respect to $W$ if $W(S)\subset S.$

Denote

$\Omega_{1}=\left\\{(x,y,z)\in\Omega:x\leq\frac{g_{0}}{g_{1}}\right\\}.$

$\Omega_{2}=\\{(x,y,z)\in\Omega_{1}:z<z^{*}\\}.$

######  Proposition 2.

The sets $\Omega_{1}$ and $\Omega_{2}$ are invariant with respect to the
operator $W$.

###### Proof.

(1) Let $(x,y,z)\in\Omega_{1}$, i.e., $g_{0}\leq x\leq\frac{g_{0}}{g_{1}}.$
Then

$x^{\prime}-\frac{g_{0}}{g_{1}}=g_{0}+x\left(1-g_{1}-cxy\right)-\frac{g_{0}}{g_{1}}\leq
g_{0}+\frac{g_{0}}{g_{1}}\left(1-g_{1}-cxy\right)-\frac{g_{0}}{g_{1}}=-\frac{g_{0}cxy}{g_{1}}\leq
0.$

Thus $(x^{\prime},y^{\prime},z^{\prime})\in\Omega_{1}$, i.e.,
$W(\Omega_{1})\subset\Omega_{1}.$

(2) Let $(x,y,z)\in\Omega_{2}$, i.e., $z<z^{*}.$ Then

$z^{\prime}=z\left(1-d_{0}+r_{1}x-r_{2}x^{2}\right)=z\left(1-r_{2}(x-\frac{r_{1}}{2r_{2}})^{2}+\frac{r_{1}^{2}-4r_{2}d_{0}}{4r_{2}}\right)\leq
z<z^{*}.$

So $(x^{\prime},y^{\prime},z^{\prime})\in\Omega_{2}$, i.e.,
$W(\Omega_{2})\subset\Omega_{2}.$

∎

Let $r_{1}^{2}=4r_{2}d_{0}$. The following theorem gives full description of
the set of limit points for the trajectory of any initial point
$\left(x^{(0)},y^{(0)},z^{(0)}\right)$ in the set $\Omega$.

######  Theorem 4.

Assume that (2.1) and (2.4) holds. For the operator $W$ given by (1.2) the
following hold:

  * (i)

If $x^{(n)}>\frac{g_{0}}{g_{1}}$ for any natural number $n$, then the
trajectory of the initial point $\left(x^{(0)},y^{(0)},z^{(0)}\right)$ taken
from the set $\Omega\setminus\Omega_{1}$, is equal to the limit point
$u^{*}_{1}$.

  * (ii)

The trajectory of any initial point $\left(x^{(0)},y^{(0)},z^{(0)}\right)$ in
the set $\Omega_{2}$ converges to fixed point $u^{*}_{1}$.

  * (iii)

If $z^{(n)}>z^{*}$ for any natural number $n$, then the trajectory of the
initial point $\left(x^{(0)},y^{(0)},z^{(0)}\right)$ in the set
$\Omega_{1}\setminus\Omega_{2}$, is equal to the limit point $u^{*}_{2}$.

###### Proof.

We have

$\begin{array}[]{lll}x^{(n)}=g_{0}+x^{(n-1)}\left(1-g_{1}-cy^{(n-1)}\right),\\\\[8.53581pt]
y^{(n)}=\frac{s_{1}\left(x^{(n-1)}\right){{}^{2}}}{s_{2}+\left(x^{(n-1)}\right){{}^{2}}}z^{(n-1)}+(1-k)y^{(n-1)},\\\\[8.53581pt]
z^{(n)}=-r_{2}\left(x^{(n-1)}-\frac{r_{1}}{2r_{2}}\right)^{2}z^{(n-1)}+z^{(n-1)}.\end{array}$
(4.4)

First, we prove the assertion (i). Let all values of $x^{(n)}$ are greater
than $\frac{g_{0}}{g_{1}}$. Then

$x^{(n)}-x^{(n-1)}=g_{0}-g_{1}x^{(n-1)}-cx^{(n-1)}y^{(n-1)}=g_{1}\left(\frac{g_{0}}{g_{1}}-x^{(n-1)}\right)-cx^{(n-1)}y^{(n-1)}\leq
0.$

$z^{(n)}-z^{(n-1)}=-r_{2}\left(x^{(n-1)}-\frac{r_{1}}{2r_{2}}\right)^{2}z^{(n-1)}\leq
0.$

Therefore, both sequences $x^{(n)}$ and $z^{(n)}$ are decreasing. Since both
sequences $x^{(n)}$ and $z^{(n)}$ are decreasing and bounded from below, we
have:

$\lim\limits_{n\to\infty}x^{(n)}\geq\frac{g_{0}}{g_{1}},\
\lim\limits_{n\to\infty}z^{(n)}\geq 0.$ (4.5)

We estimate $x^{(n)}$ and $z^{(n)}$ by the following:

$x^{(n)}=g_{0}-g_{1}x^{(n-1)}-cx^{(n-1)}y^{(n-1)}+x^{(n-1)}<g_{0}+(1-g_{1})x^{(n-1)}$

$<g_{0}+(1-g_{1})\left(g_{0}+(1-g_{1})x^{(n-2)}\right)$

$<g_{0}+g_{0}(1-g_{1})+...+g_{0}(1-g_{1})^{n-1}+(1-g_{1})^{n}x^{(0)}$

$=\frac{g_{0}}{g_{1}}\left(1-(1-g_{1})^{n}\right)+(1-g_{1})^{n}x^{(0)}.$

$z^{(n)}=\left(1-r_{2}\left(x^{(n-1)}-\frac{r_{1}}{2r_{2}}\right)^{2}\right)z^{(n-1)}<\left(1-r_{2}\left(\frac{g_{0}}{g_{1}}-\frac{r_{1}}{2r_{2}}\right)^{2}\right)z^{(n-1)}$

$<\left(1-r_{2}\left(\frac{g_{0}}{g_{1}}-\frac{r_{1}}{2r_{2}}\right)^{2}\right)^{2}z^{(n-2)}<...<\left(1-r_{2}\left(\frac{g_{0}}{g_{1}}-\frac{r_{1}}{2r_{2}}\right)^{2}\right)^{n}z^{(0)}.$

Thus
$x^{(n)}<\frac{g_{0}}{g_{1}}\left(1-(1-g_{1})^{n}\right)+(1-g_{1})^{n}x^{(0)},\
z^{(n)}<\left(1-r_{2}\left(\frac{g_{0}}{g_{1}}-\frac{r_{1}}{2r_{2}}\right)^{2}\right)^{n}z^{(0)}.$
Consequently

$\lim\limits_{n\to\infty}x^{(n)}\leq\frac{g_{0}}{g_{1}},\
\lim\limits_{n\to\infty}z^{(n)}\leq 0.$ (4.6)

Based on the inequalities (4.5) and (4.6), we can conclude that the sequence
$x^{(n)}$ converges to $\frac{g_{0}}{g_{1}}$, and the sequence $z^{(n)}$
converges to 0, respectively. From (4.4) it follows
$\lim\limits_{n\to\infty}y^{(n)}=0.$

Let’s prove the assertion (ii). Let
$\left(x^{(0)},y^{(0)},z^{(0)}\right)\in\Omega_{2}.$ Since $\Omega_{2}$ is an
invariant set, then $g_{0}\leq x^{(n)}\leq\frac{g_{0}}{g_{1}}$, $0\leq
y^{(n)}\leq C$, $0\leq z^{(n)}<z^{*}$ for any natural number $n.$ Moreover,
since $z^{(n)}$ is strictly decreasing,
$0\leq\lim\limits_{n\to\infty}z^{(n)}<z^{*}.$ Suppose
$\lim\limits_{n\to\infty}z^{(n)}=\overline{z}\neq 0$. From the third equation
of system (4.4) it is clear that the existence of the limit $z^{(n)}$ implies
the existence of the limit $x^{(n)}$. Similarly, from the first equation of
system (4.4), the existence of the limit $x^{(n)}$ implies the existence of
the limit $y^{(n)}$. We denote the limits of the sequences $x^{(n)}$ and
$y^{(n)}$ by $\overline{x}$ and $\overline{y}$, respectively. Then from the
system (4.4), we form the following:

$\left\\{\begin{array}[]{lll}\overline{x}=g_{0}+\overline{x}(1-g_{1}-c\overline{y}),\\\\[2.84526pt]
\overline{y}=\frac{s_{1}\overline{x}{{}^{2}}}{s_{2}+\overline{x}{{}^{2}}}\overline{z}+(1-k)\overline{y},\\\\[2.84526pt]
\overline{z}=\left(1-r_{2}\left(\overline{x}-\frac{r_{1}}{2r_{2}}\right)^{2}\right)\overline{z}.\end{array}\right.$
(4.7)

By (4.7) we have

$\overline{x}=\frac{r_{1}}{2r_{2}}=x^{*},\ \overline{y}=y^{*},\
\overline{z}=z^{*}$ (see (2.2)).

From this contradiction it follows that the sequence $z^{(n)}$ converges to 0.
Due to the convergence of $z^{(n)}$ to zero, Lemma 3 (second part) guarantees
that $y^{(n)}$ also converges to zero. Leveraging this result and Lemma 3
(third part), we can further conclude that $x^{(n)}$ converges to
$g_{0}/g_{1}$.

(iii). Let all values of $z^{(n)}$ are greater than $z^{*}$. Then, since
$z^{(n)}$ is decreasing and bounded from below,
$\lim\limits_{n\to\infty}z^{(n)}\geq z^{*}.$ From the third equation of system
(4.4) it is clear that the existence of the limit $z^{(n)}$ implies the
existence of the limit $x^{(n)}$. Similarly, from the first equation of system
(4.4), the existence of the limit $x^{(n)}$ implies the existence of the limit
$y^{(n)}$. By (4.4) we obtain

$\lim\limits_{n\to\infty}x^{(n)}=x^{*},\
\lim\limits_{n\to\infty}y^{(n)}=y^{*},\
\lim\limits_{n\to\infty}z^{(n)}=z^{*}.$

∎

Figure 1. Considering the system defined by equation (1.2) with the following
parameter values:
$g_{0}=0.1,g_{1}=0.15,c=0.6,s_{1}=1,s_{2}=0.2,k=0.1,r_{1}=r_{2}=1$, and
$d_{0}=0.4$, where $r_{1}^{2}<4r_{2}d_{0}$ holds. We observe that for both
initial points, (0.5, 0.3, 0.016) and (0.5, 0.18, 0.016), the trajectory of
the system converges to the fixed point (2/3, 0, 0). Figure 2.  Considering
the system defined by equation (1.2) with the following parameter values:
$g_{0}=0.1,g_{1}=0.15,c=0.6,s_{1}=1,s_{2}=0.2,k=0.1,r_{1}=r_{2}=1$, and
$d_{0}=0.25$, where $r_{1}^{2}=4r_{2}d_{0}$ holds.We observe the following
behavior: For the initial point (0.5, 0.3, 0.016): The trajectory of the
system converges to the fixed point (2/3, 0, 0). For the initial point (0.5,
0.18, 0.016): The trajectory of the system converges to the fixed point (1/2,
1/12, 3/200).

###  4.1. Biological interpretation

(1.2) model represents the dynamics of glucose, insulin, and $\beta$-cell
mass. In this section, we predict the normal behavior of the glucose
regulatory system and the pathways leading to diabetes, depending on the
parameter values and initial conditions.

Each point (vector) $(G,I,\beta)=(x,y,z)\in\Omega$ can be considered as a
state (a measure) of glucose, insulin, and $\beta$-cell mass. It can be seen
from Proposition 1 that system (1.2) has two fixed points $u^{*}_{1}$ and
$u^{*}_{2}$. The fixed point $u^{*}_{1}$ corresponds to the disease pathology
in the model, while $u^{*}_{2}$ represents a physiological state with healthy
glucose levels. $d_{0}$ is the death rate at zero glucose.

Let us give some interpretations of our main results:

  * •

Assume that (2.1) holds and let $r_{1}^{2}<4r_{2}d_{0}$. Under this condition
on $d_{0}$ (i.e. the death rate at zero glucose), the value of $\beta$-cell
mass and the level of insulin decrease, and the level of glucose increases. In
summary, regardless of the initial state $(G,I,\beta)$ of the system (1.2), it
converges to a fixed point, $u^{*}_{1}$, which represents a pathological
condition (Interpretation of Theorem 3).

  * •

(Case Theorem 4, part i.) Assume that (2.1) holds and let
$r_{1}^{2}=4r_{2}d_{0}$. In addition, even if the glucose level is always
greater than $g_{0}/g_{1}$, the system approaches a pathological state.

  * •

(Case Theorem 4, part ii.)If the value of $\beta$ mass in the initial state is
less than $z^{*}$, then the system approaches the pathological state.

  * •

(Case Theorem 4, part iii.) If the value of $\beta$ mass is always greater
than $z^{*}$, then the system approaches the physiological state.

## References

  * [1] Lyubich Y.I. (1992), Mathematical structures in population genetics, Springer-Verlag, Berlin.
  * [2] Topp B., Promislow K., Devries G., Miura R. M., Finegood D. T (2000), A model of $\beta$-cell mass, insulin, and glucose kinetics: pathways to diabetes, Journal of Theoretical Biology 206 (4) 60-619.
  * [3] Devaney R.L. (2003), An Introduction to Chaotic Dynamical System (Westview Press).
  * [4] Wang Ch, Li X., (2014), Stability and Neimark-Sacker bifurcation of a semi-discrete population model, Journal of Applied Analysis and Computation, vol.4, no.4, 419-435
  * [5] Goel P., Insulin resistance or hypersecretion? the $\beta IG$ picture revisited, Journal of Theoretical Biology 384 (2015) 131-139. https://doi.org/10.1016/j.jtbi.2015.07.033 PMID: 26300065
  * [6] Rozikov U.A. (2020), Population dynamics: algebraic and probabilistic approach. World Sci. Publ. Singapore, 460 pp.
  * [7] Boxonov Z.S., Rozikov U.A. (2021), A discrete-time dynamical system of stage-structured wild and sterile mosquito population, Nonlinear studies, Vol. 28, No.2, p.413-425.
  * [8] Boxonov Z.S., Rozikov U.A. (2021), Dynamical system of a mosquito population with distinct birth-death rates, Journal of Appleid Nonlinear Dynamics, Vol.10, No.4, p.807-816.
  * [9] Boxonov Z.S. (2023),A discrete-time dynamical system of mosquito population, JDEA., Vol. 29, No. 1, p.67-83. doi.org/10.1080/10236198.2022.2160245
  * [10] Yang B.,, Li J., Haller M.J., Schatz D.A., Rong L. (2023), Modeling the progression of Type 2 diabetes with underlying obesity, https://doi.org/10.1371/journal.pcbi.1010914


# Model Predictive Control based Energy Management System for Home Energy
Resiliency

1st Ninad Gaikwad  School of Electrical Engineering & Computer Science
Washington State University
Pullman, USA
<EMAIL_ADDRESS>2nd Shishir Lamichhane  School of Electrical
Engineering & Computer Science
Washington State University
Pullman, USA
<EMAIL_ADDRESS>3rd Anamika Dubey  School of Electrical
Engineering & Computer Science
Washington State University
Pullman, USA
<EMAIL_ADDRESS>Ninad Gaikwad, Shishir Lamichhane, and Anamika Dubey
School of Electrical Engineering & Computer Science
Washington State University
Pullman, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>and<EMAIL_ADDRESS>
###### Abstract

As the occurrence of extreme weather events is increasing so are the outages
caused by them. During such unplanned outages, a house needs to be provided
with an energy supply to maintain habitable conditions by maintaining thermal
comfort and servicing at least critical loads. An energy system consisting of
rooftop photovoltaic (PV) panels along with battery storage is an excellent
carbon-free choice to provide energy resiliency to houses against extreme
weather-related outages. However, to provide habitable conditions this energy
system has to provide not only for the non-air-conditioning (non-AC) load
demand but also for the turning on of the AC system which has a considerably
higher startup power requirement as compared to its rated power. Hence, an
intelligent automated decision-making controller is needed which can manage
the trade-off between competing requirements.

In this paper, we propose such an intelligent controller based on Model
Predictive Control (MPC). We compare its performance with a Baseline
controller which is unintelligent, and a Rule-Based controller which has some
intelligence, based on three resiliency metrics that we have developed. We
perform extensive simulations for numerous scenarios involving different
energy system sizes and AC startup power requirements. Every simulation is one
week long and is carried out for a single-family detached house located in
Florida in the aftermath of Hurricane Irma in 2017. The simulation results
show that the MPC controller performs better than the other controllers in the
more energy-constrained scenarios (smaller PV-battery size, larger AC startup
power requirement) in providing both thermal comfort and servicing non-AC
loads in a balanced manner.

###### Index Terms:

Energy Resiliency against Extreme Weather Events, Home Energy Management
System, Off-grid Home Energy Management, Model Predictive Control, PV Battery
Systems, Home Air Conditioning

##  I Introduction

There has been a substantial increase in the occurrence of extreme weather
events in recent years, where the average number of such events causing
damages exceeding $ 1 billion has increased from 8.1 each year to 22 events
per year for the last three years (2021-2023) [1]. Hurricane Idalia (2023),
the great Texas freeze (2021), California flooding (2023), Southern heatwaves
(2023), and Hawaii wildfires (2023) are some well-known extreme weather events
that occurred in the recent past that caused immense damages. In addition,
such events can greatly damage the power systems infrastructure thereby
causing long-unplanned power outages which lead to diverse socio-economic
damages [2]. For example, the power system network of Puerto Rico and the US
Virgin Islands experienced significant damages in 2017 as a result of
Hurricanes Irma and Maria which led to the longest blackout, lasting almost 11
months leaving 3.3 million Puerto-Ricans without power [3].

A house equipped with rooftop solar PV together with battery storage systems
(HPVB) to provide an excellent carbon-free alternative to provide electricity
supply to residential customers during such an unplanned outage to maintain
habitable conditions until the grid is restored [4]. Holistic habitable
conditions can be maintained by providing some level of thermal comfort along
with servicing the critical loads (refrigerators, freezers, cooking ranges,
etc.) almost all the time, and servicing the non-critical loads whenever
possible. Such exigencies for maintaining habitable conditions lead to a need
for a Home Energy Management System (HEMS) that can manage the trade-off
between competing requirements optimally.

A plethora of work has been done to develop HEMS for the grid-connected case
(no outage) where the primary objective is to increase energy savings and/or
decrease energy costs, see [5], [6], [7], and [8], but this is not the focus
of our paper. On the contrary, the amount of work done to develop HEMS that
can provide energy resiliency against extreme weather events is limited and
can be broadly divided into two: planning-based approaches and operation-based
approaches. The planning-based approaches in [9], and [10] plan for optimal
sizing of the PV and battery system so that the dual objectives of reducing
the cost of the energy system and providing resiliency against outages are
met. However, as these approaches do not provide any intelligent operational
support during the outage there can arise situations where the planned system
size is inadequate to provide for habitable conditions owing to a multitude of
reasons like longer than expected outages, lower than expected solar
irradiance, and changes in the expected house demand patterns.

Now, to reliably provide energy resiliency to a house, operation-based
approaches which provide explicit control decisions to optimally manage the
home energy system during outages, have been developed to mitigate the
shortcomings of the planning-based approaches. In our previous works [11], and
[12] we have developed operation-based HEMS that are based on MPC and
Reinforcement Learning (RL) respectively for an HPVB; but, we have not
considered AC operation, a multi-level load prioritization scheme, or real
house load demand data. Moreover, [13] has developed an MPC-based HEMS for an
HPVB similar to our work [11]; but, it does not consider AC operation, any
discrete on-off decisions, or a realistic load prioritization scheme. Also,
[14] has developed a rule-based HEMS for an HPVB with an additional Electric
Vehicle (EV) battery system; but, it does not consider AC operation, circuit-
level load grouping, real house load demand data, or the effect of variation
of PV and battery size. Finally [15] is the closest to our work to the best of
our knowledge, it has developed an HEMS based on MPC for an HPVB; but, it does
not consider the higher startup power requirement of the AC, a multi-level
load prioritization scheme, any discrete on-off decisions, the effect of
variation of PV and battery size, or a multi-dimensional resiliency metric to
evaluate the controller performance.

In this work, we develop an operation-based HEMS based on an MPC architecture
to provide holistic home energy resilience against long-unplanned power
outages caused by extreme weather events. The major contributions of this
work: 1. We consider the operation of AC along with its higher startup power
requirement (more realistic problem formulation), 2. We consider a flexible
multi-level load prioritization scheme that is based on circuit-level load
grouping (easier real-world implementation, see the products and
services/products of SPAN.io, Inc.), 3. We consider the discrete AC on-off
decisions explicitly in our optimization formulation for the MPC, while the
on-off decisions for the circuit-level prioritized loads are computed by a
separate priority stack controller (making the resulting mixed integer linear
program computationally simpler), 4. We assess the performance of our proposed
controller against two other controllers based on a three-dimensional
resiliency metric, 5. We consider the effect of different PV, and battery
system sizes and AC startup power requirements while assessing the performance
of our proposed controller, and 6. We use real house load demands, and weather
data after a real extreme weather event for our simulations.

The rest of the paper is organized as follows: the system is described in
Section II, the plant model that the controllers try to manage is explained in
Section III, the mathematical formulations for the proposed MPC controller
along with the baseline controller and rule-based smart controller are laid
out in Section IV the simulation setup used to generate results is discussed
in Section V, the simulation results are presented and discussed in Section
VI, and the main conclusions are provided in Section VII.

##  II System Description

Figure 1: Schematic of the home energy system with its load description.

The Figure 1 illustrates the home energy system with a detailed classification
of its intrinsic loads (energy-consuming systems). The house is equipped with
a solar PV system for energy generation and a battery system for energy
storage. The two broad categories of the loads in the house are the AC load
and non-AC loads referred to as other loads. The other loads are further
categorized based on their priority level which depends on their subjective
importance for maintaining a habitable condition for a given house based on
the owner’s preferences during a grid outage (we consider that grid outage is
caused by an extreme weather event and it is not known when the grid will be
restored). It should be noted that the priority levels and the loads that are
part of it as shown in Figure 1 are not absolute but a preference of the
authors. The objective of any well-designed controller for this system will be
to serve as much load as possible to maintain maximum comfort conditions with
the constraint of finite energy by controlling the charging/discharging of the
battery system, the on-off of the AC system, and the on-off of the prioritized
other loads. Furthermore, two sensors are present to measure the internal
house temperature and the battery energy level so that controllers can use
these measurements as feedback if and when required. Also, we consider that
the on-off of the AC system and that of the prioritized other loads are
realized through a smart switchboard (each circuit corresponding to both the
AC load and individual other loads has a smart switch associated with them in
the smart switchboard).

##  III Plant Model

The plant model is a mathematical description of the real system, using which
a digital simulator can be developed for testing the performance of a
controller (being designed for managing the real system) in a closed-loop
simulation. Now for describing our system (presented in Section II)
mathematically we consider time to be discrete, $k=0,1,2,3,\dots,$ denotes the
time index, and $\Delta T_{s}$ is the time between any two adjacent time
indices (i.e time-step). Moreover, $E(k)$ is any energy produced or consumed
by any subsystem between the time interval marked by time indices $k-1$ and
$k$. The dependence on $k$ for any variable can be omitted where necessary and
constants will be mentioned explicitly. Any variables with a bar on top
($\bar{.}$) or bar on the bottom ($\underaccent{\bar}{.}$) are either
constants or data (i.e. not computed but provided). The units for the upcoming
energy, power, and temperature variables are $kWh$, $kW$, and ${}^{\circ}C$
respectively.

In our system, the maximum energy available from the PV system
($\bar{E}_{pv}$) is modeled using a set of equations given in [11] and [16],
and it is affected by weather variables: solar irradiance
($\bar{G}HI\;(kW/m^{2})$), ambient outside temperature ($\bar{T}_{am}$), and
wind speed ($\bar{W}_{s}\;(m/s)$). In addition, it has two dynamic systems:
the house thermal dynamics system that gives the time evolution of $T_{h}$
(internal house temperature), and the battery energy dynamics system that
gives the time evolution of $E_{bat}$ (battery energy level). These are
modeled as linear discrete-time dynamical systems and are presented in detail
in [17] and [11] respectively. The control commands for our system are given
as follows: battery system control command is $u_{bat}=[c,d]^{T}$ where
$c,d\in\\{0,1\\}$ (binary variables) are the battery charging and discharging
commands respectively, AC system on-off control command is
$u_{ac}\in\\{0,1\\}$ (binary variable), and other loads on-off control
commands are given as $u_{l_{i}}\in\\{0,1\\}\;\forall i\in I=\\{1,\dots,n\\}$
(binary variables) where $n$ is the number of different prioritized loads and
the set $I$ is ordered in descending order of priority. The energy consumed by
AC now is given as $u_{ac}\bar{E}_{ac}$ where $\bar{E}_{ac}$ is rated energy
consumption of the AC system, and that consumed by the other loads is given by
$E_{l}=\sum_{i=1}^{n}u_{l_{i}}\bar{E}_{l_{i}}$ and we have
$\bar{E}_{l}=\sum_{i=1}^{n}\bar{E}_{l_{i}}$ where $\bar{E}_{l_{i}}$ is the
energy demanded by the $i^{\text{th}}$ priority load and $\bar{E}_{l}$ is
total energy demand from the other loads ($i\in I$). The remaining plant model
consists of equations representing the energy consumption/production
interaction between the Battery system, PV system, AC system, and other loads.
These equations are similar to the ones described in [11] with an additional
power mismatch constraint over the energy mismatch constraint (both are
explained in Section IV-C).

##  IV Control Algorithms and Performance Metrics

(a) MPC Controller Schematic.

(b) Rule-Based Smart Controller Schematic.

Figure 2: Schematics of MPC and Rule-Based Smart Controllers..

###  IV-A Model Predictive Controller

The MPC illustrated in Figure 2(a) computes the control commands at discrete
time indices $k=1,2\dots N$ where $N$ is the total number of time indices in
the planning horizon (only the first time index i.e the $j^{\text{th}}$ time
index control commands are sent to the plant), for the next time index $j+1$,
the computation is done taking into account the feedback of the dynamic states
of the plant and the new available forecasts/estimates of the disturbances).
The decision variables for the optimization problem elemental to the MPC
controller are as follows: the states of the process
$x(k)=[T_{h}(k),E_{bat}(k)]^{T}$; the control commands
$u(k)=[\Gamma(k),u_{ac}(k),E_{l}(k)]^{T}$, where $\Gamma(k)$ (converted to
realizable control commands as shown in (1) and (2)), $u_{ac}(k)$ and
$E_{l}(k)$ (converted to realizable control commands by the Priority Stack
Controller described in Section IV-A4) are the fraction of the normal battery
charging energy, the AC on-off control command and energy consumed by the
other loads respectively (it should be noted that having
$\Gamma(k)\in\mathbb{R}$ eliminates 2 binary variables, while having
$E_{l}(k)$ eliminates $n$ binary variables); the internal variables
$v(k)=[E_{pv}(k),\zeta_{h}(k),\zeta_{l}(k),f_{on}(k),f_{off}(k),\theta_{bat}(k)]^{T}$,
where $E_{pv}(k)$ is the energy produced by the PV panels, $\zeta_{h}(k)$ and
$\zeta_{l}(k)$ are the slack variable for house temperature and minimum
critical load ($\bar{E}_{cri}\triangleq\bar{E}_{l_{1}}$ i.e. the highest
priority load) respectively to ensure feasibility, $f_{on}(k)$ and
$f_{off}(k)$ are the AC turn-on and turn-off indicators as described in Eq.(3)
and Eq.(4), and $\theta_{bat}(k)$ is the battery discharge indicator. The
exogenous inputs whose predictions are assumed to be known for the $N$ time
step $w(k)=[\bar{E}_{pv}(k),\bar{E}_{l}(k),\bar{T}_{am}(k)]^{T}$, where
$\bar{E}_{pv}$, $\bar{E}_{l}(k)$, and $\bar{T}_{am}$ are the available energy
from the PV panels, the forecast of the total other load demand of the house,
and the ambient outside air temperature respectively. Hence, the complete
decision vector for the optimization problem is given as $[X,U,V]^{T}$, where
$X:=[x(k+1),\dots,x(k+N)]^{T}$, $U:=[u(k),\dots,u(k+N-1)]^{T}$ and
$V:=[v(k),\dots,v(k+N-1)]^{T}$. Moreover; $u_{ac}$, $f_{on}$, $f_{off}$ and
$\theta_{bat}$ are binary variables i.e. $\in{0,1}$ while rest all variables
$\in\mathbb{R}$.

$\displaystyle c(k)$ $\displaystyle=\begin{cases}1\;,&\text{if
}\Gamma(k)>0,\\\ 0\;,&\text{if }\Gamma(k)\leq 0,\\\ \end{cases}$ (1)
$\displaystyle d(k)$ $\displaystyle=\begin{cases}1\;,&\text{if
}\Gamma(k)<0,\\\ 0\;,&\text{if }\Gamma(k)\geq 0,\\\ \end{cases}$ (2)
$\displaystyle f_{on}(k)$
$\displaystyle=\begin{cases}1\;,&u_{ac}(k)=1\;\&\;u_{ac}(k-1)=0,\\\
0\;,&\text{otherwise,}\end{cases}$ (3) $\displaystyle f_{off}(k)$
$\displaystyle=\begin{cases}1\;,&u_{ac}(k)=0\;\&\;u_{ac}(k-1)=1,\\\
0\;,&\text{otherwise.}\end{cases}$ (4)

####  IV-A1 Optimization Problem Formulation

The MPC solves a Mixed Integer Linear Programming (MILP) problem to compute
the control commands over the planning horizon, and an instance of this
problem at any time index $j$ is given mathematically as follows:

$\displaystyle\min_{X,U,V}$
$\displaystyle\sum\limits^{j+N-1}_{k=j}\bigg{[}\lambda_{1}(N-k)\zeta_{h}(k)+\lambda_{2}(N-k)\zeta_{l}(k)-$
(5) $\displaystyle\quad\quad\lambda_{3}(N-k)E_{l}(k)-\lambda_{4}E_{bat}(k)+$
$\displaystyle\quad\quad\lambda_{5}\theta_{bat}(k)\bigg{]},$

subject to the following constraints:

$\displaystyle T_{h}(k+1)=AT_{h}(k)+Bu_{ac}(k)Q_{ac}+D\bar{T}_{am}(k),$
$\displaystyle E_{bat}(k+1)=E_{bat}(k)-\Gamma(k)\bar{E}_{bat}^{c,dc},$
$\displaystyle
u_{ac}(k)\bar{E}_{ac}-\Gamma(k)\bar{E}_{bat}^{c,dc}+E_{l}(k)=E_{pv}(k),$ (8)
$\displaystyle
u_{ac}(k)=\sum\limits^{k}_{m=j}f_{on}(m)-\sum\limits^{k}_{m=j}f_{off}(m),$ (9)
$\displaystyle
f_{on}(k)\bar{P}_{ac}\leq\theta_{bat}(k)\bar{P}_{bat}^{dc}+\frac{\bar{E}_{pv}(k)}{\Delta
T_{s}},$ (10) $\displaystyle\underaccent{\bar}{T}_{h}\leq
T_{h}(k)\leq\bar{T}_{h}+\zeta_{h}(k),$ (11)
$\displaystyle\underaccent{\bar}{E}_{bat}\leq E_{bat}(k)\leq\bar{E}_{bat},$
(12) $\displaystyle\underaccent{\bar}{\Gamma}\leq\Gamma(k)\leq\bar{\Gamma},$
(13) $\displaystyle\bar{E}_{cri}(k)-\zeta_{l}(k)\leq
E_{l}(k)\leq\bar{E}_{l}(k),$ (14) $\displaystyle 0\leq
E_{pv}(k)\leq\bar{E}_{pv}(k),$ (15) $\displaystyle f_{on}(k)+f_{off}(k)\leq
1,$ (16) $\displaystyle\zeta_{h}(k)\geq 0,$ (17)
$\displaystyle\theta_{bat}(k)\geq\Gamma(k),$ (18)
$\displaystyle\bar{E}_{cri}(k)\geq\zeta_{l}(k)\geq 0.$ (19)

####  IV-A2 Cost Function Description

The cost function is given in (5), where the first term in the cost function
tries to minimize the temperature slack variable $\zeta_{h}$, which implicitly
tries to keep the house temperature within bounds, see (11). The second term
tries to minimize energy demand slack variable $\zeta_{l}$, which implicitly
tries to maximize the energy provision for critical load demand, see (14). The
third term tries to maximize the house energy demand variable $E_{h}$
(excluding AC energy demand) as we would like most of the house energy demand
to be met most of the time. The fourth term tries to maximize the battery
energy level $E_{bat}$, this is to ensure that the battery is charged when
there is excess solar and that it is maintained at an appropriate level to
provide for future house energy demand. The final term tries to minimize the
battery charging mode indicator $\theta_{bat}$, this minimization coupled with
the constraint in (18) helps carve out the right behavior of this variable
i.e. it takes the value of 1 when the battery is in charging mode and take the
value 0 when the battery is in discharging mode. The $\lambda$s are weighting
factors for the cost function terms to achieve the right trade-off between the
competing objectives resulting in desired control behavior. The $(N-k)$ terms
present with certain cost function terms are also additional weighting factors
to incorporate temporal weighting which descends through the planning horizon.
The particular objective with this temporal weighting is desired to be
achieved for the present and near future but is deemed okay to fail in the far
future w.r.t. the planning horizon, as this is a finite energy system we
desire certain of our objectives to be achieved myopically/greedily w.r.t.
time.

####  IV-A3 Constraints Description

(LABEL:eq:House_eqcon_1) and (LABEL:eq:Battery_eqcon_1) are the house thermal
dynamics and the battery dynamics respectively (there is no plant-model
mismatch for house thermal dynamics while there is a plant-model mismatch for
battery dynamics as we have modeled $\Gamma$ instead of $c$ and $d$ and
ignored the charging/discharging efficiencies). (8) is the energy balance
equation (there is a plant-model mismatch for energy balance as we have
modeled $E_{l}$ instead of $\sum_{i=1}^{n}u_{l_{i}}\bar{E}_{l_{i}}$). (9) is
the relationship of the AC on-off control variable with the binary indicators
for the turning on ($f_{on}$) and turning off ($f_{off}$) of the AC. (10)
conveys the AC startup power constraint where $f_{on}$ indicates the turning
on of AC while the binary indicator $\theta_{bat}$ indicates the charging mode
of the battery. (11) gives the bounds (upper bound $\bar{T}_{h}$ and lower
bound $\underaccent{\bar}{T}_{h}$) for the house temperature where $\zeta_{h}$
serves as the slack variable for the upper bound of the house temperature (we
require it due to finite energy supply of the system). (12) gives the bounds
(upper bound $\bar{E}_{bat}$ and lower bound $\underaccent{\bar}{E}_{bat}$) on
the battery energy storage, while (13) gives the bounds (upper bound
$\bar{\Gamma}$ and lower bound $\underaccent{\bar}{\Gamma}$) on the
charging/discharging rate of the battery. (14) provides upper and lower bounds
for the load energy demand variable where $E_{cri}$ is the energy demanded by
the critical loads, $\bar{E}_{l}$ is the actual total energy demand of the
house (AC excluded), and $\zeta_{l}$ is the slack variable (we require it due
to finite energy supply of the system).(15) gives the bound on the PV energy
production where $\bar{E}_{pv}$ is the maximum energy production potential of
the installed PV. (16) shows the complementary nature of the AC turning on
($f_{on}$) and turning off ($f_{off}$) binary indicators. (17) and (19)
provide for the slack variables $\zeta_{h}$ and $\zeta_{l}$ to be positive
respectively. (18) pits the battery charging mode indicator $\theta_{bat}$
against the battery charging/discharging rate so that it takes the value of 1
only when the battery is in the charging mode (we have added minimization of
$\theta_{bat}$ to the cost function for appropriate behavior).

####  IV-A4 Priority Stack Controller

The priority stack controller converts the first time index ($j$ w.r.t. the
MILP planning horizon) control command value $E_{l}(j)$ into realizable
control commands $u_{l_{m}}(j)\;,m\in I$ which can be sent to the plant. The
priority stack controller logic is given as follows;

$\displaystyle u_{l_{m}}(j)$ $\displaystyle=\begin{cases}1\;,&\text{if
}\sum_{i=1}^{m}\bar{E}_{l_{m}}(j)\leq E_{l}(j),\\\ 0\;,&\text{otherwise,}\\\
\end{cases}$ (20)

and the complete control command for the priority stack controller can be
given as $u_{l}^{ps}=[u_{l_{1}},\dots,u_{l_{m}},\dots,u_{l_{n}}]^{T}$.

###  IV-B Baseline Controller

Figure 3: Schematic of the Baseline Controller.

The baseline controller illustrated in Figure 3 is meant to abstract the
operation of a commercially installed PV-Battery system for a house in the
absence of any intelligent/supervisory control during a grid outage. It
consists of three controllers: the battery control logic, the AC thermostat
controller, and the load control logic. It should be noted these three
controllers operate independently of each other at the level of their
concerned subsystem, and that a smart switchboard is not present.

####  IV-B1 Battery Control Logic

The battery control logic commands the charging/discharging of the battery and
is given as follows:

$\displaystyle c(k)=\begin{cases}1,&\text{if }\bar{E}_{pv}(k)\geq E_{d}(k)\\\
0,&\text{otherwise }\\\ \end{cases}$ (21) $\displaystyle
d(k)=\begin{cases}1,&\text{if }E_{pv}(k)<E_{d}(k)\\\ 0,&\text{otherwise,}\\\
\end{cases}$ (22)

where $E_{d}$ is the total house energy demand given by
$E_{d}(k)=E_{l}(k)+u_{ac}(k)\bar{E}_{ac}$. This controller is a mathematical
abstraction of the interaction of the battery system with the rest of the
subsystems in the house along with the logic present in the charge controller
which connects the battery to the inverter (the battery charging/discharging
is not controlled by an external controller but rather is a manifestation from
the intrinsic logic and circuit conditions).

####  IV-B2 AC Thermostat Controller

The AC thermostat controller commands the on-off of the AC and is given as
follows:

$\displaystyle u_{ac}(k)=\begin{cases}1,&\text{if }T_{h}(k)\geq\bar{T}_{h}\\\
0,&\text{if }T_{h}(k)\leq\underaccent{\bar}{T}_{h}\\\
u_{ac}(k-1),&\text{otherwise.}\end{cases}$ (23)

This is a realized controller that is intrinsic to the AC system.

####  IV-B3 Load Control Logic

The load logic controller controls the on-off of the other loads and is given
as follows:

$\displaystyle u_{l_{i}}(k)=\begin{cases}1,&\text{if }\bar{E}_{l_{i}}(k)>0\\\
0,&\text{if }\bar{E}_{l_{i}}(k)=0.\end{cases}$ (24)

This controller is a mathematical abstraction of the physical turning on and
turning off of the individual other loads by any occupant of the house.

###  IV-C Rule-Based Controller

Algorithm 1 Algorithm for the Rule-Based Controller

Use the Solar PV model to compute $\bar{E}_{pv}(k)$

Use the Baseline Controller (21)-(24) to compute $u^{b}(k)$.

Compute $E_{bat_{d}}(k)$ as given in (25).

Compute $S_{ac}^{on}(k)$.

Compute $E_{mis}(k)$ using (28).

Compute $P_{mis}(k)$ using (31).

if ($E_{mis}(k)\geq 0$) AND ($P_{mis}(k)\geq 0$) then

$u(k)\leftarrow u^{b}(k)$

else if ($(E_{mis}(k)\geq 0$) AND ($P_{mis}(k)<0$) then

$u_{ac}(k)\leftarrow 0$

$u_{bat}(k)\leftarrow u^{b}_{bat}(k)$

$u_{l}(k)\leftarrow u^{b}_{l}(k)$

else if $(E_{mis}(k)<0)$ then

if $(|E_{mis}(k)|\leq u^{b}_{ac}(k)\bar{E}_{ac})$ then

$u_{ac}(k)\leftarrow 0$

$u_{bat}(k)\leftarrow u^{b}_{bat}(k)$

$u_{l}(k)\leftarrow u^{b}_{l}(k)$

else if $(|E_{mis}(k)|\leq u^{b}_{ac}(k)\bar{E}_{ac})+\bar{E}_{l}(k)$ then

$u_{ac}(k)\leftarrow 0$

$u_{bat}(k)\leftarrow u^{b}_{bat}(k)$

Compute $u^{ps}_{l}(k)$ using the Priority Stack Controller given in (20).

$u_{l}(k)\leftarrow u^{ps}_{l}(k)$

end if

end if

The rule-based controller illustrated in Figure 2(b) is developed to check if
a mathematically simpler controller based on just sophisticated rules can
provide a performance close to the mathematically dense MPC controller.
Algorithm 1 presents the algorithm for the rule-based controller at an
arbitrary time index $k$.

This algorithm returns the complete control command
$u(k)=[u_{bat}(k),u_{ac}(k),u_{l}(k)]^{T}$, $u_{bat}(k)=[c(k),d(k)]^{T}$ is
the battery control command, and
$u_{l}(k)=[u_{l_{1}}(k),\dots,u_{l_{n}}(k)]^{T}$ is the other load command.
Similarly, $u^{b}(k)=[u^{b}_{bat}(k),u^{b}_{ac}(k),u^{b}_{l}(k)]^{T}$ is the
complete control command computed by the Baseline Controller (the Baseline
Controller logic is an intrinsic part of this controller),
$u^{b}_{bat}(k)=[c^{b}(k),d^{b}(k)]^{T}$ is the battery control command, and
$u^{b}_{l}(k)=[u^{b}_{l_{1}}(k),\dots,u^{b}_{l_{n}}(k)]^{T}$ is the other load
command.

The battery charging/discharging dispatch ($E_{bat_{c}}(k)$ and
$E_{bat_{d}}(k)$ respectively) are computed as follows;

$\displaystyle E_{bat_{d}}(k)$
$\displaystyle=d^{b}(k)min(E_{bat}^{c,dc},\bar{E}_{bat}(k)-E_{bat}).$ (25)

$S_{ac}^{on}(k)$ is a binary indicator for the turning on of the AC (it should
be noted that this requires the controller to have a memory of the past on-off
state of the AC and this is not shown in Figure 2(b) for maintaining
simplicity). We compute the energy mismatch of the system ($E_{mis}(k)$) as
follows;

$\displaystyle E_{d}^{b}(k)$
$\displaystyle=u_{ac}^{b}(k)\bar{E}_{ac}+\bar{E}_{l}(k),$ (26) $\displaystyle
E_{a}^{b}(k)$ $\displaystyle=\bar{E}_{pv}(k)+E_{bat_{d}}(k),$ (27)
$\displaystyle\bar{E}_{mis}(k)$ $\displaystyle=E_{a}^{b}(k)-E_{d}^{b}(k),$
(28)

where $E_{d}^{b}(k)$ and $E_{a}^{b}(k)$ are the total house load demand and
total available house generation under the action of the baseline controller.
Finally, the power mismatch $P_{mis}(k)$ is computed as follows;

$\displaystyle P_{d}^{b}(k)$ $\displaystyle=S_{ac}^{on}(k)\bar{P}_{ac},$ (29)
$\displaystyle P_{a}^{b}(k)$ $\displaystyle=\dfrac{\bar{E}_{pv}(k)}{\Delta
T_{s}}+\bar{P}_{bat}(k),$ (30) $\displaystyle\bar{P}_{mis}(k)$
$\displaystyle=P_{a}^{b}(k)-\bar{P}_{d}^{b}(k),$ (31)

where $P_{d}^{b}(k)$ (we do not consider the power demand contribution of the
other loads here as the AC startup power is relatively very large, and the AC
startup transient time is very small) and $P_{a}^{b}(k)$ are the AC startup
power demanded and total available startup power under the action of the
baseline controller. It is important to note that both the Baseline and the
Rule-Based controllers are reactive as they do not plan in the presence of
forecast information.

##  V Simulation Setup

We perform closed-loop simulations of the plant model along with each of the
three controllers described in Section LABEL:section:ControlAlgorithms to
compare their performance and generate results. The setup for these
simulations is described as follows:

###  V-A Simulation Period and Location

The simulation is carried out for seven days with a time-step of 10 minutes
($\Delta T_{s}=10$ minutes or $\Delta T_{s}=1/6$ hours). The period of
simulation selected is from Sept. 11, 2017, to Sept. 17, 2017, and the
location of the house being simulated is selected to be in Gainesville, FL,
USA as hurricane Irma made landfall there on Sept. 11, 2017, and a grid outage
was experienced by a large section of its population for a few days.

###  V-B Weather and Load Demand Data

We have obtained the weather data (i.e. $\bar{GHI}$, $\bar{W}_{s}$, and
$\bar{T}_{am}$) for the concerned location and simulation period from the
National Solar Radiation Database, see [18]. The other load demand data (i.e.
$\bar{E}_{l_{i}},\;\forall i\in I$ excluding thermal load data) is obtained
from Pecan Street Dataport, see [19]. We use the load data of a house located
in Austin, TX, USA that has both PV and battery installed from the freely
available part (has a collection of 75 houses located in Austin in Texas
state, California, and New York) of the Pecan Street dataset as it is the
closest location to the simulation location (assuming the selected house will
be experiencing similar weather to the simulation location). The eight
prioritized load classifications as shown in Figure 1 stem from the various
load demands recorded and available in the Pecan Street dataset.

###  V-C Simulation Parameters

The PV system parameters are as follows: a PV module (Tesla SC325) with a
rated power of $0.325\;kW$ is used to form a $\alpha_{pv}\times 10.075\;kW$
rooftop PV system, where $\alpha_{pv}\in A_{pv}=\\{0.25,0.5,0.75,1\\}$ is the
PV size factor. The Battery system parameters are as follows: a single Tesla
Powerwall is used as the battery system with $\bar{E}_{bat}=\alpha_{bat}\times
13.5\;kWh$, $\underaccent{\bar}{E}_{bat}=0\;kWh$ (we allow for complete
utilization of storage capacity due to a finite energy system),
$\bar{P}_{bat}^{c,dc}=\alpha_{bat}\times 5\;kW$ is the rated battery
charging/discharging rate, where
$\bar{E}_{bat}^{c,dc}=\bar{P}_{bat}^{c,dc}\times\Delta T_{s}$, and
$\bar{P}_{bat}^{dc}=\alpha_{bat}\times 7\;kW$ is a short time (assumed to be
in the order of AC startup transient time) maximum battery discharging rate,
where $\alpha_{bat}\in A_{bat}=\\{0.25,0.5,0.75,1\\}$ is the Battery size
factor. It should be noted that we haven’t sized the system according to the
selected house’s load but used the most commercially viable PV-Battery system
size ($\alpha_{pv}=1$ and $\alpha_{bat}=1$) for a single-family detached house
and then progressively reduced the system size to see the effect of resource
constraint on the performance of the three controllers.

The AC system parameters are as follows: $\bar{P}_{ac}^{rated}=3\;kW$ is the
rated power of the AC, where $\bar{E}_{ac}=\bar{P}_{ac}^{rated}\times\Delta
T_{s}$, $\bar{T}_{h}=25\;^{\circ}C$, $\underaccent{\bar}{T}_{h}=23\;^{\circ}C$
(the reason for wide dead-band is to conserve energy in a finite energy
system), the AC startup power is given as
$\bar{P}_{ac}=(1-\alpha_{V})\times\alpha_{I}\times P_{ac}^{rated}$ where
$\alpha_{V}=0.3$ is the AC startup voltage factor which is almost a constant
for different AC systems and $\alpha_{I}\in A_{I}=\\{3,4,5,6,7,8\\}$ is the AC
startup current factor (i.e. locked-rotor current factor) which is variable
w.r.t. the type of AC system, see [20]. Finally, the house thermal system
parameters are used from [17].

The MPC parameters are as follows: $N=144$ i.e the planning horizon is of 24
hours with a time-step of 10 minutes, $\lambda_{i}=1\;\forall
i\in\\{1,2,3,4,5\\}$, $\bar{\Gamma}=1$, and $\underaccent{\bar}{\Gamma}=-1$.

###  V-D Computation

The computing system used for performing the simulations is a Windows-based
Desktop machine with 32Gb RAM and 3.7 GHz $\times$ 10 CPU. The plant along
with the three control algorithms are simulated in MATLAB. The MILP described
in (5)-(19) for the planning horizon of 24 hours ($N=144$) leads to an
optimization problem with $576$ binary variables and $1440$ continuous
variables.

The optimization solver used to solve the MILP formulated for the MPC
controller is GUROBI [21]. The GUROBI solver parameters that are changed from
their defaults are as follows: MIPGap = 1% (the default value of 0.01% leads
to solver stalling, hence a more pragmatic value is needed given the real-time
control requirement), MIPFocus = 2 (the default value is 0 which utilizes a
higher solution strategy to strike a balance between finding feasible
solutions and proving their optimality, we select the value 2 that focuses on
proving optimality as we need the most optimal solution given a finite energy
system), Cuts = 2 (the default value is -1 which automatically selects how the
cuts should be generated, we select the value 2 that generates cuts
aggressively for a more tightly bound solution), Presolve = 2 (the default
value is -1 which automatically selects the presolve strategy, we select the
value 2 that aggressively presolves for getting a tighter model before branch
and bound can start given the real-time control requirements), and TimeLimit =
500 seconds (the default is no time limit, the value of 500 seconds is smaller
than the simulation time-step of 10 minutes i.e. 600 seconds which helps in
the case of a simulation iteration where the solver is stalling and we can
still have some suboptimal solution at the end of 500 seconds to provide to
the plant and continue on to the next iteration). With these parameters, the
solver didn’t stall for any of the MPC closed-loop simulations and the average
time to solve the MILP was $0.29\;s$.

##  VI Results And Discussion

###  VI-A Performance Metrics

####  VI-A1 Critical Load Resiliency Metric $(LRM_{cri})$

It measures the amount of critical loads served and is the ratio of total
critical loads served to the total critical loads demanded by the system.

$\displaystyle LRM_{cri}$
$\displaystyle=\frac{\sum\limits^{T_{sim}}_{k=1}E_{cri}(k)}{\sum\limits^{T_{sim}}_{k=1}\bar{E}_{cri}(k)},$
(32)

where $T_{sim}$ is the total number of time steps in a given simulation.

####  VI-A2 Other Loads Resiliency Metric $(LRM_{o})$

It measures the amount of other loads served and is given by the ratio of
other loads served to the total amount of other loads over the entire time.

$\displaystyle LRM_{o}$
$\displaystyle=\frac{\sum\limits^{T_{sim}}_{k=1}E_{l}(k)}{\sum\limits^{T_{sim}}_{k=1}\bar{E}_{l}(k)},$
(33)

####  VI-A3 Thermal Resiliency Metric $(TRM_{h})$

It measures the average thermal comfort and is given by the ratio of the
number of times the house temperature is above its upper bound to the total
number of simulation time steps.

$\displaystyle TRM_{h}$
$\displaystyle=\frac{\sum_{k=1}^{T_{sim}}\left[\bar{T}_{h}-T_{h}(k)\right]_{+}}{T_{sim}},$
(34)

###  VI-B Comparison of Controllers for all Simulation cases $\in A_{I}\times
A_{pv}\times A_{bat}$

(a) $LRM_{cri}$ Metric.

(b) $LRM_{o}$ Metric.

(c) $TRM_{h}$ Metric.

Figure 4: Performance comparison of the three controllers based on three
resiliency metrics, for all the considered simulation cases $\in A_{I}\times
A_{pv}\times A_{bat}$.

The Fig. 4 illustrates the performance comparison of the three controllers
(i.e. MPC, Baseline, and Rule-Based) based on the three metrics (i.e.
$LRM_{cri}$, $LRM_{o}$, and $TRM_{h}$) for all the considered combinations of
the PV size, Battery size, and AC startup current factor (i.e. simulation
scenarios).

The Fig. 4(a) shows the variation of critical load resiliency metric
($LRM_{cri}$) for the three different controllers across all simulation
scenarios. It can be observed across all simulation scenarios, that the MPC
and Rule-Based controllers perform almost equivalently, and better than the
Baseline controller in providing for the critical loads as they are more
sophisticated as compared to the Baseline controller. Moreover, the
performance of all three controllers improves with an increase in the PV and
Battery size and it decreases with an increase in the AC startup factor as
expected (excluding the simulation scenarios where the performance metric
values saturate due to physical limitations on both the best and worst
performance that can be extracted from the system). On the other hand, for a
given PV size the performance of the three controllers does not improve
significantly for an increase in the Battery size, as an appropriate size of
PV is required to generate enough power during the daytime to charge the
battery to its capacity. Also, it can be observed that for smaller sizes of PV
and Battery the MPC controller performs better than the Rule-Based controller
in most scenarios, as its utilization of planning and forecasts within its
formulation becomes more important for better performance in more energy-
constrained scenarios when compared to only reactive control provided by the
Rule-Based controller.

The Fig. 4(b) shows the variation of other load resiliency metric ($LRM_{o}$)
for the three different controllers across all simulation scenarios. The
observations that can be made here are equivalent to those made for the
critical load resiliency metric, except that the relative numerical values of
$LRM_{cri}$ are greater than that of $LRM_{o}$, as critical loads are one of
the eight prioritized load components which make up the other load and have
the highest priority in the case of MPC and Rule-Based controllers.

Fig. 4(c) shows the variation of thermal resiliency metric ($TRM_{h}$) for the
three controllers across all simulation scenarios. It can be observed that the
MPC controller can maintain thermal comfort better than the Baseline and Rule-
Based controllers in most of the simulation scenarios. This is due to the
presence of a competing objective to maintain thermal comfort in the MPC
controller’s cost function which lends to some relative priority (can be
controlled by setting the right $\lambda_{1}$) given to the operation of AC
for providing thermal comfort; while no such priority for AC operation is
provided in the Rule-Based and Baseline controllers leading to their poor
performance. Other observations are equivalent to those made for the critical
($LRM_{cri}$) and other load ($LRM_{o}$) resiliency metrics, with the
exception that for some simulation cases the Baseline controller performs
better than the MPC and Rule-Based controllers in providing thermal comfort
(one such simulation scenario is explored in detail in the next subsection).
Also, it should be noted that the $\lambda_{i}$s of the MPC are not tuned for
each of the simulation scenarios individually lending to a more conservative
performance.

###  VI-C Comparison of Controllers for the Simulation case:
$\alpha_{I}=4,\;\alpha_{pv}=0.5,\;\alpha_{bat}=0.5$

Fig. 5(a), Fig. 5(b) and Fig. 5(c) illustrate the time series plots of the
desired and serviced loads due to MPC, Baseline, and Rule-Based controllers
respectively for the entire simulation period. It can be observed that all
three controllers can provide for at most the desired critical load as the
size of the PV/Battery system is not adequate to provide for all of the
desired load; except for the initial simulation period where they provide for
the desired load greater than the critical load as the house temperature is
within bounds (turning on of the AC is not required) and the battery state-of-
charge (SoC) is high (enough stored energy). Table I shows that the MPC and
Rule-Based controllers outperform the Baseline controller in providing for the
desired loads, this is expected as the Baseline controller is not as
sophisticated in load management as the other two controllers. Furthermore,
the MPC and Rule-Based controllers perform equivalently in providing for the
other loads, but the MPC controller is slightly worse in providing for the
critical loads as it prioritizes the thermal load to some extent as opposed to
the Rule-Based controller which completely sheds the thermal load in favor of
critical/other loads in an energy-constrained situation.

TABLE I: Performance metrics for the simulation case: $\alpha_{I}=4,\alpha_{pv}=0.5,\alpha_{bat}=0.5.$ | $LRM_{cri}$ | $LRM_{o}$ | $TRM_{h}$
---|---|---|---
MPC | 0.65 | 0.53 | 0.46
Baseline | 0.45 | 0.37 | 0.55
Rule-Based | 0.66 | 0.53 | 0.36

(a) MPC.

(b) Baseline.

(c) Rule-Based.

(d) MPC.

(e) Baseline.

(f) Rule-Based.

(g) MPC.

(h) Baseline.

(i) Rule-Based.

Figure 5: Performance comparison of the three controllers for serving the
other load along with the critical load (top row), thermal load (middle row),
and PV energy utilization (bottom row) for the simulation case:
$\alpha_{I}=4,\alpha_{pv}=0.5,\alpha_{bat}=0.5.$.

Fig. 5(d), Fig. 5(e) and Fig. 5(f) illustrate the time series plots of the
house temperature and battery SoC due to MPC, Baseline, and Rule-Based
controllers respectively for the entire simulation period. From these time
series plots along with the metric values presented in Table I it can be
observed that the Baseline controller outperforms the other two controllers in
providing thermal comfort, but due to lack of any energy management logic, it
prioritizes the thermal and other loads equally and in this simulation
scenario the thermal load turns out as the winner. On the other hand, the MPC
controller performs significantly better than the Rule-Based controller in
providing thermal comfort and maintaining battery SoC as it has a mathematical
formulation to prioritize both, while the Rule-Based controller lacks such a
formulation leading to a dead battery (approximately between 72 to 150 hours)
which leads to no power availability for AC startup (still the critical/other
loads are being served). Around the 72, 96, 120, and 144-hour marks where the
temperature is out of bound for both the MPC and Baseline controllers, it can
be observed that in the case of MPC, it is a planning-based decision in the
presence of forecast information to not serve the thermal load and instead
prioritize the charging of the battery and serving of the other load; while in
the case of Baseline controller the thermal load along with the other load is
not served as trying to turn on the AC causes a power limit violation leading
to no power transfer in the system. Fig. 5(g), Fig. 5(h) and Fig. 5(i)
illustrate the time series plots of the PV energy utilization against the
available PV energy for MPC, Baseline, and Rule-Based controllers respectively
for the entire simulation period, these time series plots corroborate the
previous observations.

Hence, from all the results and the associated discussion, it is clear that
the MPC controller with its capability for intelligent planning in the
presence of forecast information is well-suited for providing energy
resiliency to a house during an unplanned outage caused by extreme weather.

##  VII Conclusion

We present an MPC-based controller for providing energy resiliency to a house
equipped with roof-top solar PV and battery energy storage during unplanned
outages caused by extreme weather events while maintaining thermal comfort
through intelligent AC operation. We compare our proposed MPC controller with
a Baseline and a Rule-Based controller based on three resiliency metrics. The
simulation results show that generally, the performance of the three
controllers is directly proportional to the PV and battery size while being
inversely proportional to the AC startup factor. The MPC controller performs
better than the other controllers in the more energy-constrained scenarios
(smaller PV-battery size, larger AC startup factor) in providing both thermal
comfort and servicing critical/other loads in a balanced manner. This study
illustrates that a controller based on intelligent planning and forecast
information utilization is more reliable than reactive rule-based controllers
in providing home energy resiliency against extreme weather events.

In the future, we plan to explore: the performance of the MPC controller
during extreme cold events, the sensitivity of the MPC controller to forecast
errors, the sizing of PV and battery based on resiliency requirements, the
development of a RL based controller to avoid solving the MILP in the MPC
controller, and the development of a distributed control framework to provide
home energy resilience against extreme weather events for a community of
heterogeneous (can have different distributed energy resources) houses.

## Acknowledgment

This work is funded by the NSF award 2208783.

## References

  * [1] “U.s. billion-dollar weather and climate disasters (2023).” [Online]. Available: https://www.ncei.noaa.gov/access/billions/
  * [2] V. Do, H. McBrien, N. Flores, A. Northrop, J. Schlegelmilch, M. Kiang, and J. Casey, “Spatiotemporal distribution of power outages with climate events and social vulnerability in the usa,” _Nature Communications_ , vol. 14, p. 2470, 2023.
  * [3] “2017 hurricane season: Federal support for electricity grid restoration in the u.s. virgin islands and puerto rico.” [Online]. Available: https://www.gao.gov/products/gao-19-296
  * [4] M. Mango, J. A. Casey, and D. Hernandez, “Resilient power: A home-based electricity generation and storage solution for the medically vulnerable during climate-induced power outages,” _Futures_ , vol. 128, p. 102707, 2021\.
  * [5] H. Shareef, M. S. Ahmed, A. Mohamed, and E. Al Hassan, “Review on home energy management system considering demand responses, smart technologies, and intelligent controllers,” _Ieee Access_ , vol. 6, pp. 24 498–24 509, 2018\.
  * [6] A. Q. Badar and A. Anvari-Moghaddam, “Smart home energy management system–a review,” _Advances in Building Energy Research_ , vol. 16, no. 1, pp. 118–143, 2022.
  * [7] R. Balakrishnan and V. Geetha, “Review on home energy management system,” _Materials Today: Proceedings_ , vol. 47, pp. 144–150, 2021.
  * [8] I. Gomes, K. Bot, M. G. Ruano, and A. Ruano, “Recent techniques used in home energy management systems: A review,” _Energies_ , vol. 15, no. 8, p. 2866, 2022.
  * [9] E. Chatterji and M. D. Bazilian, “Battery storage for resilient homes,” _IEEE Access_ , vol. 8, pp. 184 497–184 511, 2020.
  * [10] E. Chatterji, K. Anderson, and M. D. Bazilian, “Planning for a resilient home electricity supply system,” _IEEE Access_ , vol. 9, pp. 133 774–133 785, 2021.
  * [11] N. Gaikwad, N. S. Raman, and P. Barooah, “Smart home energy management system for power system resiliency,” in _2020 IEEE Conference on Control Technology and Applications (CCTA)_. IEEE, 2020, pp. 1072–1079.
  * [12] N. S. Raman, N. Gaikwad, P. Barooah, and S. P. Meyn, “Reinforcement learning-based home energy management system for resiliency,” in _2021 American Control Conference (ACC)_. IEEE, 2021, pp. 1358–1364.
  * [13] P. Haessig, R. Bourdais, H. Guéguen _et al._ , “Resilience in energy management system: A study case,” _IFAC-PapersOnLine_ , vol. 52, no. 4, pp. 395–400, 2019.
  * [14] A. K. Candan, A. R. Boynuegri, and N. Onat, “Home energy management system for enhancing grid resiliency in post-disaster recovery period using electric vehicle,” _Sustainable Energy, Grids and Networks_ , vol. 34, p. 101015, 2023\.
  * [15] M. Wei, W. Jia, and Y. Fu, “Energy management and control system for a pv-battery system to improve residential building resiliency under extreme weather conditions.”
  * [16] D. Faiman, “Assessing the outdoor operating temperature of photovoltaic modules,” _Progress in Photovoltaics: Research and Applications_ , vol. 16, no. 4, pp. 307–315, 2008.
  * [17] B. Cui, C. Fan, J. Munk, N. Mao, F. Xiao, J. Dong, and T. Kuruganti, “A hybrid building thermal modeling approach for predicting temperatures in typical, detached, two-story houses,” _Applied energy_ , vol. 236, pp. 101–116, 2019\.
  * [18] National Renewable Energy Laboratory (NREL), “Nsrdb: National solar radiation database,” https://nsrdb.nrel.gov/data-viewer, accessed: 2024-02-04.
  * [19] Pecan Street Inc., “Pecan street dataport,” https://www.pecanstreet.org/dataport/, accessed: 2024-02-04.
  * [20] Lazar Rozenblat, “Locked-rotor amps (lra) calculator,” https://generators.smps.us/surge-current-calculator.html, accessed: 2024-02-04.
  * [21] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2023. [Online]. Available: https://www.gurobi.com