text:
\begin{document}

\title[$P-$extremal]{\bf A global domination principle for $P-$pluripotential theory}

\author{Norm Levenberg* and Menuja Perera}{\thanks{*Supported by Simons Foundation grant No. 354549}}
\subjclass{32U15, \ 32U20, \ 31C15}
\keywords{$P-$extremal function, global domination principle}

\maketitle

\begin{center} {\bf In honor of 60 years of Tom Ransford} \end{center}
\begin{abstract} We prove a global domination principle in the setting of $P-$pluripotential theory. This has many applications including a general product property for $P-$extremal functions. The key ingredient is the proof of the existence of a strictly plurisubharmonic $P-$potential.
\end{abstract}

\section{Introduction}\label{sec:intro} Following \cite{Bay}, in \cite{BBL} and \cite{BosLev} a pluripotential theory associated to plurisubharmonic (psh) functions on $\CC^d$ having growth at infinity specified by $H_P(z):=\phi_P(\log |z_1|,...,\log |z_d|)$ where
$$\phi_P(x_1,...,x_d):=\sup_{(y_1,...,y_d)\in P}(x_1y_1+\cdots + x_dy_d)$$
is the indicator function of a convex body $P\subset (\RR^+)^d$ was developed. Given $P$, the classes
$$L_P=L_P(\CC^d):= \{u\in PSH(\CC^d): u(z)- H_P(z) =0(1), \ |z| \to \infty \}$$
and
$$L_P^+ =L_P^+(\CC^d):=\{u\in L_P: u(z)\geq H_P(z) +c_u\}$$
are of fundamental importance. These are generalizations of the standard Lelong classes
$L(\CC^d)$, the set of all plurisubharmonic (psh) functions $u$ on $\CC^d$ with $u(z) - \max[\log |z_1|,...,\log |z_d|] = 0(1), \ |z| \to \infty$, and
$$L^+(\CC^d)=\{u\in L(\CC^d): u(z)\geq \max[0,\log |z_1|,...,\log |z_d|]  + C_u\}$$
which correspond to $P=\Sigma$ where
$$\Sigma:=\{(x_1,...,x_d)\in \RR^d: x_1,...,x_d \geq 0, \ x_1+\cdots + x_d\leq 1\}.$$
For more on standard pluripotential theory, cf., \cite{[K]}.

Given $E\subset \CC^d$, the {\it $P-$extremal function of $E$} is defined as $V^*_{P,E}(z):=\limsup_{\zeta \to z}V_{P,E}(\zeta)$ where
$$V_{P,E}(z):=\sup \{u(z):u\in L_P(\CC^d), \ u\leq 0 \ \hbox{on} \ E\}.$$
For $P=\Sigma$, we write $V_E:=V_{\Sigma,E}$. For $E$ bounded and nonpluripolar, $V_E^*\in L^+(\CC^d)$; $V_E^*=0$ q.e. on $E$ (i.e., on all of $E$ except perhaps a pluripolar set); and $(dd^c V_E^*)^d=0$ outside of $\bar E$ where $(dd^c V_E^*)^d$ is the complex Monge-Amp\`ere measure of $V_E^*$ (see section 2). A key ingredient in verifying a candidate function $v\in L^+(\CC^d)$ is equal to $V_E^*$ is the following global domination principle of Bedford and Taylor:

\begin{proposition} \label{domprinc2} \cite{BT88} Let $u\in L(\CC^d)$ and $v\in L^+(\CC^d)$ and suppose $u \leq v$ a.e.-$(dd^c v)^d$. Then $u\leq v$ on $\CC^d$.
\end{proposition}

Thus if one finds $v\in L^+(\CC^d)$ with $v=0$ a.e. on $\bar E$ and $(dd^cv)^d=0$ outside of $\bar E$ then $v=V_E^*$. For the proof of Proposition \ref{domprinc2} in \cite{BT88} the fact that in the definition of the Lelong classes $\max[\log |z_1|,...,\log |z_d|] $ and $\max[0,\log |z_1|,...,\log |z_d|] $ can be replaced by the K\"ahler potential
$$u_0(z):=\frac{1}{2}\log {(1+|z|^2)}:= \frac{1}{2}\log (1+\sum_{j=1}^d |z_j|^2)$$ is crucial; this latter function is strictly psh and $(dd^cu_0)^d>0$ on $\CC^d$.

We prove a version of the global domination principle for very general $L_P$ and $L_P^+$ classes. We consider convex bodies $P\subset (\RR^+)^d$ satisfying
\begin{equation}\label{phyp} \Sigma \subset kP \ \hbox{for some} \ k\in \ZZ^+.\end{equation}

\begin{proposition} \label{gctp} For $P\subset (\RR^+)^d$ satisfying (\ref{phyp}), let $u \in L_P$ and $v \in L_P^+$ with $u \leq v$ a.e. $MA(v)$. Then $u \leq v$ in $\CC^d$.
\end{proposition}

As a corollary, we obtain a generalization of Proposition 2.4 of \cite{BosLev} on $P-$extremal functions:

\begin{proposition}  \label{productproperty} Given $P\subset (\RR^+)^d$ satisfying (\ref{phyp}), let $E_1,...,E_d\subset \CC$ be compact and nonpolar. Then
\begin{equation}\label{prodprop} V^*_{P,E_1\times \cdots \times E_d}(z_1,...,z_d)=\phi_P(V^*_{E_1}(z_1),...,V^*_{E_d}(z_d)).\end{equation}
\end{proposition}

The main issue in proving Proposition 1.2 (Proposition \ref{gctp} below) is the construction of a strictly psh $P-$potential $u_P$ which can replace the logarithmic indicator function $H_P(z)$ used to define $L_P$ and $L_P^+$. To do this, we utilize a classical result on subharmonic functions in the complex plane which we learned in Tom Ransford's beautiful book \cite{Rans}; thus it is fitting that this article is written in his honor.

\section{The global $P-$domination principle}\label{sec:back}
 Following \cite{BBL} and \cite{BosLev}, we fix a convex body $P\subset (\RR^+)^d$; i.e., a compact, convex set in $(\RR^+)^d$ with non-empty interior $P^o$. The most important example is the case where $P$ is the convex hull of a finite subset of $(\ZZ^+)^d$ in $(\RR^+)^d$ with $P^o\not =\emptyset$ ($P$ is a non-degenerate convex polytope). Another interesting class consists of the $(\RR^+)^d$ portion of an $\ell^q$ ball for $1\leq q\leq \infty$; see (\ref{ellq}). Recall that $H_P(z):=\phi_P(\log |z_1|,...,\log |z_d|)$ where $\phi_P$ is the indicator function of $P$.

A $C^2-$function $u$ on $D\subset \CC^d$ is {\it strictly psh} on $D$ if the complex Hessian $H(u):=\sum_{j,k=1}^d {\partial^2u\over \partial z_j  \partial \bar z_k}$ is positive definite on $D$. We define
$$dd^cu :=2i \sum_{j,k=1}^d {\partial^2u\over \partial z_j  \partial \bar z_k} dz_j \wedge d\bar z_k $$
and
$$(dd^cu)^d=dd^cu \wedge \cdots \wedge dd^cu=c_d \det H(u) dV$$
where $dV=({i\over 2})^d \sum_{j=1}^d dz_j\wedge d\bar z_j$ is the volume form on $\CC^d$ and $c_d$ is a dimensional constant. Thus $u$ strictly psh on $D$ implies that $(dd^cu)^d=fdV$ on $D$ where $f>0$. We remark that if $u$ is a locally bounded psh function then $(dd^cu)^d$ is well-defined as a positive measure, the {\it complex Monge-Amp\`ere measure} of $u$; this is the case, e.g., for functions $u\in L_P^+$.

\begin{definition} \label{ppot} We say that $u_P$ is a {\it strictly psh $P-$potential} if
\begin{enumerate}
\item $u_P\in L_P^+$ is strictly psh on $\CC^d$ and
\item there exists a constant $C$ such that $|u_P(z)-H_P(z)|\leq C$ for all $z\in \CC^d$.
\end{enumerate}
\end{definition}

This property implies that $u_P$ can replace $H_P$ in defining the $L_P$ and $L_P^+$ classes:
$$L_P= \{u\in PSH(\CC^d): u(z)- u_P(z) =0(1), \ |z| \to \infty \}$$
and
$$L_P^+ =\{u\in L_P: u(z)\geq u_P(z) +c_u\}.$$
Given the existence of a strictly psh $P-$potential, we can follow the proof of Proposition \ref{domprinc2} in \cite{BT88} to prove:

\begin{proposition} \label{gctp} For $P\subset (\RR^+)^d$ satisfying (\ref{phyp}), let $u \in L_P$ and $v \in L_P^+$ with $u \leq v$ a.e. $(dd^cv)^d$. Then $u \leq v$ in $\CC^d$.
\end{proposition}

\begin{proof} Suppose the result is false; i.e., there exists $z_0\in \CC^d$ with $u(z_0)>v(z_0)$. Since $v\in L_P^+$, by adding a constant to $u,v$ we may assume $v(z)\geq u_P(z)$ in $\CC^d$. Note that $(dd^c u_P)^d>0$ on $\CC^d$. Fix $\delta, \ \epsilon>0$ with $\delta <\epsilon /2$ in such a way that the set
$$S:=\{z\in \CC: u(z)+\delta u_P(z)\}>(1+\epsilon)v(z)\}$$
contains $z_0$. Then $S$ has positive Lebesgue measure. Moreover, since $\delta <\epsilon$ and $v\geq u_P$, $S$ is bounded. By the comparison principle (cf., Theorem 3.7.1 \cite{[K]}), we conclude that
$$\int_S (dd^c[u+\delta u_P])^d\leq \int_S (dd^c (1+\epsilon)v)^d.$$
But $\int_S (dd^c \delta u_P)^d>0$ since $S$ has positive Lebesgue measure, so
$$(1+\epsilon) \int_S (dd^c v)^d>0.$$
By hypothesis, for a.e.-$(dd^cv)^d$ points in supp$(dd^cv)^d\cap S$ (which is not empty since $\int_S (dd^c v)^d>0$), we have
$$(1+\epsilon)v(z)< u(z)+\delta u_P(z)\leq v(z)+\delta u_P(z),$$
i.e., $v(z)< \frac{1}{2}u_P(z)$ since $\delta <\epsilon /2$. This contradicts the normalization $v(z)\geq u_P(z)$ in $\CC^d$.
\end{proof}

In the next section, we show how to construct $u_P$ in Definition \ref{ppot} for a convex body in $(\RR^+)^d$ satisfying (\ref{phyp}).

\section{Existence of strictly psh $P-$potential}
For the $P$ we consider, $\phi_P\geq 0$ on $(\RR^+)^d$ with $\phi_P(0)=0$. We write $z^J=z_1^{j_1}\cdots z_d^{j_d}$ where  $J=(j_1,...,j_d)\in P$ (the components $j_k$ need not be integers) so that $$H_P(z):=\sup_{J\in P} \log |z^J|:=\phi_P(\log^+ |z_1|,...,\log^+ |z_d|)$$
with $|z^J|:=|z_1|^{j_1}\cdots |z_d|^{j_d}$. To construct a strictly psh $P-$potential $u_P$, we first assume $P$ is a convex polytope in $(\RR^+)^d$ satisfying (\ref{phyp}). Thus $(a_1,0,...,0), \ldots,  (0,...,0,a_d)\in \partial P$ for some $a_1,...,a_d>0$. A calculation shows that
$$\log (1 + |z_1|^{2a_1} +\cdots + |z_d|^{2a_d})$$
is strictly psh in $\CC^d$.

We claim then that
\begin{equation}\label{up} u_P(z):= \frac{1}{2} \log( 1+ \sum_{J\in Extr(P)}|z^J|^2)\end{equation}
is strictly psh in $\CC^d$ and the $L_P, \ L_P^+$ classes can be defined using $u_P$ instead of $H_P$; i.e., $u_P$ satisfies (1) and (2) of Definition \ref{ppot}. Here, $Extr(P)$ denotes the extreme points of $P$ but we omit the origin ${\bf 0}$. Note that $(a_1,0,...,0), \ldots,  (0,...,0,a_d)\in Extr(P)$.

Indeed, in this case,
$$H_P(z)=\sup_{J\in P} \log |z^J|=\max[0,\max_{J\in Extr(P)} \log |z^J|]$$
so clearly for $|z|$ large, $|u_P(z)-H_P(z)|=0(1)$ while on any compact set $K$,
$$\sup_{z\in K} |u_P(z)-H_P(z)|\leq C=C(K)$$
which gives (2) (and therefore that $u_P\in L_P^+$).

It remains to verify the strict psh of $u_P$ in (\ref{up}). We use reasoning based on a classical univariate result which is exercise 4 in section 2.6 of \cite{Rans}: if $u, v$ are nonnegative with $\log u$ and $\log v$ subharmonic (shm) -- hence $u, v$ are shm -- then
$\log (u + v)$ is shm. The usual proof is to show $(u + v)^a$ is shm for any $a > 0$ -- which is exercise 3 in section 2.6 of \cite{Rans} -- which trivially
follows since $u, v$ are shm and $a > 0$.
However, assume $u, v$ are smooth and compute the Laplacian $\Delta \log (u + v)$ on $\{u,v>0\}$:
$$\bigl(\log (u + v)\bigr)_{z\bar z} = \frac{(u+v)(u_{z\bar z} + v_{z\bar z})-(u_z+v_z)(u_{\bar z} + v_{\bar z})}{(u+v)^2}$$
$$=\frac{[uu_{z\bar z}- |u_z|^2 +vv_{z\bar z}- |v_z|^2]+[uv_{z\bar z}+vu_{z\bar z}- 2\Re (u_zv_{\bar z})]}{(u+v)^2}.$$
Now $\log u, \ \log v$ shm implies $uu_{z\bar z}- |u_z|^2\geq 0$ and $vv_{z\bar z}- |v_z|^2\geq 0$ with strict inequality in case of strict shm. Since $\log (u + v)$ is shm, the entire numerator is nonnegative:
$$[uu_{z\bar z}- |u_z|^2 +vv_{z\bar z}- |v_z|^2]+[uv_{z\bar z}+vu_{z\bar z}- 2\Re (u_zv_{\bar z})]\geq 0$$
so that the ``extra term''
$$uv_{z\bar z}+vu_{z\bar z}- 2\Re (u_zv_{\bar z})$$
is nonnegative whenever
$$(\log u)_{z\bar z} + (\log v)_{z\bar z}=uu_{z\bar z}- |u_z|^2 +vv_{z\bar z}- |v_z|^2=0.$$
We show $\Delta \log (u+v)$ is strictly positive on $\{u,v>0\}$ if one of $\log u$ or $\log v$ is strictly shm.

\begin{proposition} Let $u, v\geq 0$ with $\log u$ and $\log v$ shm. If one of $\log u$ or $\log v$ is strictly shm, e.g., $\Delta \log u > 0$,
then $\Delta \log (u+v) > 0$ on $\{u,v>0\}$.
\end{proposition}

\begin{proof} We have $u,v\geq 0$, $u_{z\bar z}, v_{z \bar z}\geq 0$, $vv_{z\bar z}- |v_z|^2\geq 0$ and $uu_{z\bar z}- |u_z|^2>0$ if $u>0$.
We want to show that
$$uv_{z\bar z}+vu_{z\bar z}- 2\Re (u_zv_{\bar z})=uv_{z\bar z}+vu_{z\bar z}-(u_zv_{\bar z}+v_zu_{\bar z})>0$$
on $\{u,v>0\}$. We start with the identity
\begin{equation}\label{tada}(uv_z-vu_z)(uv_{\bar z}-v u_{\bar z})=u^2v_zv_{\bar z}+v^2 u_z u_{\bar z}-uv(u_zv_{\bar z}+v_zu_{\bar z})\geq 0.\end{equation}
Since $uu_{z\bar z}- |u_z|^2>0$ and $vv_{z\bar z}- |v_z|^2\geq 0$,
$$uu_{z\bar z}>u_zu_{\bar z}, \ vv_{z\bar z}\geq v_z v_{\bar z}.$$
Thus
$$uv_{z\bar z}+vu_{z\bar z}=\frac{u}{v}vv_{z\bar z}+\frac{v}{u}uu_{z\bar z}>\frac{u}{v}v_zv_{\bar z}+\frac{v}{u}u_zu_{\bar z}.$$
Thus it suffices to show
$$\frac{u}{v}v_zv_{\bar z}+\frac{v}{u}u_zu_{\bar z} \geq u_zv_{\bar z}+u_{\bar z}v_z.$$
Multiplying both sides by $uv$, this becomes
$$u^2v_zv_{\bar z}+v^2u_zu_{\bar z} \geq uv(u_zv_{\bar z}+u_{\bar z}v_z).$$
This is (\ref{tada}).

\end{proof}

\noindent This proof actually shows that
$$uv_{z\bar z}+vu_{z\bar z}- 2\Re (u_zv_{\bar z})>0$$
under the hypotheses of the proposition.

\begin{remark} To be precise, this shows strict shm only on $\{u,v>0\}$. In the multivariate case, this shows the restriction of $\log (u + v)$ to the intersection of a complex line and $\{u,v>0\}$ is strictly shm if one of
$\log u, \log v$ is strictly psh so that $\log (u + v)$ is strictly psh on $\{u,v>0\}$. \end{remark}

Now with $u_P$ in (\ref{up}) we may write
$$u_P(z)=\log (u+v)$$
where
\begin{equation}\label{ueqn}u(z)=1 + |z_1|^{2a_1} +\cdots + |z_d|^{2a_d}\end{equation}
 -- so that $\log u$ is strictly psh in $\CC^d$ -- and
$$v(z)=  \sum_{J\in Extr(P)}|z^J|^2- |z_1|^{2a_1} -\cdots - |z_d|^{2a_d}.$$ If $v\equiv 0$ (e.g., if $P=\Sigma$) we are done. Otherwise $v\geq 0$ (being a sum of nonnegative terms) and $\log v$ is psh (being the logarithm of a sum of moduli squared of holomorphic functions) showing that $u_P(z):= \frac{1}{2} \log( 1+ \sum_{J\in Extr(P)}|z^J|^2)$ is strictly psh where $v>0$. There remains an issue at points where $v=0$ (coordinate axes). However, if we simply replace the decomposition $u_P(z)=\log (u+v)$ by $u_P(z)=\log (u_{\epsilon}+v_{\epsilon})$ where
$$u_{\epsilon}:= 1 + (1-\epsilon)(|z^a_1|^2+....+|z^a_d|^2) \ \hbox{and} $$
$$v_{\epsilon}:= \sum_{J\in Extr P}|z^J|^2-(1-\epsilon)(|z^a_1|^2+....+|z^a_d|^2)$$
for $\epsilon>0$ sufficiently small, then the result holds everywhere. We thank F. Piazzon for this last observation.

If $P\subset (\RR^+)^d$ is a convex body satisfying (\ref{phyp}), we can approximate $P$ by a monotone decreasing sequence of convex polytopes $P_n$ satisfying the same property. Since $P_{n+1}\subset P_n$ and $\cap_n P_n =P$, the sequence $\{u_{P_n}\}$ decreases to a function $u\in L_P^+$. Since each $u_{P_n}$ is of the form
$$u_{P_n}(z)=\log (u_n+v_n)$$
where $u_n(z)=1 + |z_1|^{2a_{n1}} +\cdots + |z_d|^{2a_{nd}}$ and $a_{nj}\geq a_j$ for all $n$ and each $j=1,...,d$ in (\ref{ueqn}), it follows that $u=:u_P$ is strictly psh and hence satisfies Definition \ref{ppot}. This concludes the proof of Proposition \ref{gctp}.

\begin{remark} Another construction of a strictly psh $P-$potential as in Definition \ref{ppot} which is based on solving a real Monge-Amp\`ere equation and which works in more general situations was recently given by C. H. Lu \cite{CL}. Indeed, his construction, combined with Corollary 3.10 of \cite{DDL}, yields a new proof of the global domination principle, Proposition \ref{gctp}.

\end{remark}

\section{The product property}

In this section, we prove the product property stated in the introduction:
\begin{proposition} \label{productproperty} For $P\subset (\RR^+)^d$ satisfying (\ref{phyp}), let $E_1,...,E_d\subset \CC$ be compact and nonpolar. Then
\begin{equation}\label{prodprop} V^*_{P,E_1\times \cdots \times E_d}(z_1,...,z_d)=\phi_P(V^*_{E_1}(z_1),...,V^*_{E_d}(z_d)).\end{equation}
\end{proposition}

\begin{remark} One can verify the formula
$$V_{P,T^d}(z)=H_P(z) =\sup_{ J\in P} \log |z^{ J}|$$
for the $P-$extremal function of the torus
$$T^d:=\{(z_1,...,z_d): |z_j|=1, \ j=1,...,d\}$$
for a general convex body by modifying the argument in \cite{[K]} for the standard extremal function of a ball in a complex norm. Indeed, let $u\in L_P$ with $u\leq 0$ on $T^d$. For $w=(w_1,...,w_d)\not \in T^d$ and $w_j\not =0$, we consider
$$v(\zeta_1,...,\zeta_d):=u(w_1/\zeta_1,...,w_d/\zeta_d) - H_P( w_1/\zeta_1,...,w_d/\zeta_d).$$
This is psh on $0< |\zeta_j|< |w_j|, \ j=1,...,d$. Since $u\in L_P$, $v$ is bounded above near the pluripolar set given by the union of the coordinate planes in this polydisk and hence extends to the full polydisk. On the boundary $|\zeta_j|= |w_j|$, $v\leq 0$ so at $(1,1,...,1)$ we get $u(w_1,...,w_d) \leq H_P( w_1,...,w_d)$. Note
$$H_P(z)=\sup_{ J\in P} \log |z^{ J}|=\phi_P(\log^+|z_1|,...,\log^+|z_d|)$$
and $V_{T^1}(\zeta)=\log^+|\zeta|$ so this is a special case of Proposition \ref{productproperty}.
\end{remark}

\begin{proof}
For simplicity we consider the case $d=2$ with variables $(z,w)$ on $\CC^2$. As in \cite{BosLev}, we may assume $V_E$ and $V_F$ are continuous. Also, by approximation we may assume $\phi_P$ is smooth. We write
$$v(z,w):=\phi_P(V_E(z),V_F(w)).$$
An important remark is that, since $P\subset (\RR^+)^2$, $P$ is convex, and $P$ contains $k\Sigma$ for some $k>0$, the function $\phi_P$ on $(\RR^+)^2$ satisfies
\begin{enumerate}
\item $\phi_P\geq 0$ and $\phi_P(x,y)=0$ only for $x=y=0$;
\item $\phi_P$ is nondecreasing in each variable; i.e., $(\phi_P)_x, (\phi_P)_y\geq 0$;
\item $\phi_P$ is convex; i.e., the real Hessian $H_{\RR}(\phi_P)$ of $\phi_P$ is positive semidefinite; and, more precisely, by the homogenity of $\phi_P$; i.e., $\phi_P(tx,ty)=t\phi_P(x,y)$,
$$\det H_{\RR}(\phi_P)=0 \ \hbox{away from the origin.}$$
\end{enumerate}

As in \cite{BosLev}, to see that
$$v(z,w)\leq V_{P,E\times F}(z,w),$$
since $\phi_P(0,0)=0$, it suffices to show that $\phi_P(V_{E}(z),V_{F}(w))\in L_P(\CC^2)$. From the definition of $\phi_P$,
$$\phi_P(V_{E}(z),V_{F}(w))=\sup_{(x,y)\in P} [xV_E(z)+yV_F(w)]$$
which is a locally bounded above upper envelope of plurisubharmonic functions. As $\phi_P$ is convex and $V_E, V_F$ are continuous, $\phi_P(V_{E}(z),V_{F}(w))$ is continuous. Since $V_E(z)=\log |z|+0(1)$ as $|z|\to \infty$ and $V_F(w)=\log |w|+0(1)$ as $|w|\to \infty$, it follows that
$\phi_P(V_{E}(z),V_{F}(w))\in L_P(\CC^2)$.

By Proposition \ref{gctp}, it remains to show $(dd^cv)^2=0$ outside of $E\times F$. Since we can approximate $v$ from above uniformly by a decreasing sequence of smooth psh functions by convolving $v$ with a smooth bump function, we assume $v$ is smooth and compute the following derivatives:
$$v_z=(\phi_P)_x (V_E)_z, \ v_w=(\phi_P)_y (V_F)_w;$$
$$v_{z\bar z} = (\phi_P)_{xx} |(V_E)_z|^2+ (\phi_P)_x (V_E)_{z\bar z};$$
$$v_{z\bar w} = (\phi_P)_{xy} (V_E)_z (V_F)_{\bar w};$$
$$v_{w\bar w} = (\phi_P)_{yy} |(V_F)_w|^2+ (\phi_P)_y (V_F)_{w\bar w}.$$
It follows from (2) that $v_{z\bar z}, v_{w\bar w}\geq 0$. Next, we compute the determinant of the complex Hessian of $v$ on $(\CC\setminus E)\times (\CC\setminus F)$ (so $(V_E)_{z\bar z}=(V_F)_{w\bar w}=0$):
$$v_{z\bar z}v_{w\bar w}-|v_{z\bar w}|^2$$
$$=(\phi_P)_{xx} |(V_E)_z|^2(\phi_P)_{yy} |(V_F)_w|^2-[(\phi_P)_{xy}]^2 |(V_E)_z|^2 |(V_F)_w|^2=$$
$$=|(V_E)_z|^2|(V_F)_w|^2 [(\phi_P)_{xx}(\phi_P)_{yy}-(\phi_P)_{xy}]^2].$$
This is nonnegative by the convexity of $\phi_P$ and, indeed, it vanishes on $(\CC\setminus E)\times (\CC\setminus F)$
by (3). The general formula for the determinant of the complex Hessian of $v$ is
$$v_{z\bar z}v_{w\bar w}-|v_{z\bar w}|^2$$
$$=|(V_E)_z|^2|(V_F)_w|^2 [(\phi_P)_{xx}(\phi_P)_{yy}-(\phi_P)_{xy}]^2]+(\phi_P)_{xx} |(V_E)_z|^2  (\phi_P)_y (V_F)_{w\bar w}$$
$$+ (\phi_P)_{yy} |(V_F)_w|^2(\phi_P)_x (V_E)_{z\bar z}+(\phi_P)_x (V_E)_{z\bar z}(\phi_P)_y (V_F)_{w\bar w}.$$
If, e.g., $z\in E$ and $w\in  (\CC\setminus F)$,
$$|(V_E)_z|^2|(V_F)_w|^2 [(\phi_P)_{xx}(\phi_P)_{yy}-(\phi_P)_{xy}]^2]=0$$
by (3) (since $(V_E(z),V_F(w))=(0,a)\not =(0,0)$) and $(V_F)_{w\bar w}=0$ so
$$v_{z\bar z}v_{w\bar w}-|v_{z\bar w}|^2=(\phi_P)_{yy} |(V_F)_w|^2(\phi_P)_x (V_E)_{z\bar z}.$$
However, we claim that
$$(\phi_P)_{yy}(0,a)=0 \ \hbox{if} \ a >0$$
since we have $\phi_P (0,ty)=t\phi_P(0,y)$. Hence
$$v_{z\bar z}v_{w\bar w}-|v_{z\bar w}|^2=0$$
if $z\in E$ and $w\in  (\CC\setminus F)$. Similarly,
$$(\phi_P)_{xx}(a,0)=0 \ \hbox{if} \ a >0$$
so that
$$v_{z\bar z}v_{w\bar w}-|v_{z\bar w}|^2=0$$
if $z\in (\CC\setminus E)$ and $w\in F$.

\end{proof}

\begin{remark} In \cite{BosLev}, a (much different) proof of Proposition \ref{productproperty} was given under the additional hypothesis that $P\subset (\RR^+)^d$ be a {\it lower set}: for each $n=1,2,...$, whenever $(j_1,...,j_d) \in nP\cap (\ZZ^+)^d$ we have $(k_1,...,k_d) \in nP\cap (\ZZ^+)^d$ for all $k_l\leq j_l, \ l=1,...,d$.

\end{remark}

Finally, although computation of the $P-$extremal function of a product set is now rather straightforward, even qualitative properties of the corresponding Monge-Amp\`ere measure are less clear. To be concrete, for $q\ge1$, let
\begin{equation}\label{ellq}
P_q:=\{(x_1,...,x_d): x_1,..., x_d\geq 0, \ x_1^q+\cdots +x_d^q\leq 1\}
\end{equation}
be the $(\RR^+)^d$ portion of an $\ell^q$ ball. Then for $1/q'+1/q=1$ we have $\phi_{P_q}(x)=||x||_{\ell^{q'}}$ (for $q =\infty$ we take $q'  =1$ and vice-versa). Hence if $E_1,...,E_d\subset \CC$,
\begin{align*}
V^*_{P_q,E_1\times \cdots \times E_d}(z_1,...,z_d)&=\|[V^*_{E_1}(z_1),V^*_{E_2}(z_2),
\ldots V^*_{E_d}(z_d)]\|_{\ell^{q'}}\\
&=
[V^*_{E_1}(z_1)^{q'}+\cdots +V^*_{E_d}(z_d)^{q'}]^{1/q'}.
\end{align*}

In the standard case $q=1$, $P_1=\Sigma$ and we have the well-known result that
$$V^*_{E_1\times \cdots \times E_d}(z_1,...,z_d)=\max[V^*_{E_1}(z_1),V^*_{E_2}(z_2),
\ldots, V^*_{E_d}(z_d)].$$
Then if none of the sets $E_j$ are polar,
$$(dd^c V^*_{E_1\times \cdots \times E_d})^d=\mu_{E_1}\times \cdots \times \mu_{E_d}$$
where $\mu_{E_j}=\Delta V^*_{E_j}$ is the classical equilibrium measure of $E_j$.

\begin{question} What can one say about supp$(dd^c V^*_{P_q,E_1\times \cdots \times E_d})^d$ in the case when $q>1$?
\end{question}

As examples, for $T^d=\{(z_1,...,z_d): |z_j|=1, \ j=1,...,d\}$ we have $V_T(z_j)=\log^+|z_j|$ and hence for $q\geq 1$
$$V_{P_q,T^d}(z)=\phi_{P_q}(\log^+|z_1|,...,\log^+|z_d|)=[\sum_{j=1}^d (\log^+|z_j|)^{q'}]^{1/q'}.$$
The measure $(dd^cV_{P_q,T^d})^d$ is easily seen to be invariant under the torus action and hence is a positive constant times Haar measure on $T^d$. Thus in this case supp$(dd^cV_{P_q,T^d})^d=T^d$ for $q\geq 1$.

For the set $[-1,1]^d$ we have $V_{[-1,1]}(z_j)=\log|z_j+\sqrt{z_j^2-1}|$ and hence for $q\geq 1$
$$
 V_{P_q,[-1,1]^d}(z_1,...,z_d)=\left\{\sum_{j=1}^d \left(\log\left|z_j+\sqrt{z_j^2-1}\right|\right)^{q'}\right\}^{1/q'}.
$$
In this case, it is not clear for $q>1$ whether supp $(dd^c  V_{P_q,[-1,1]^d})^d=[-1,1]^d$.


text:
\begin{document}

\maketitle

\begin{abstract}

In this article, the focus is mainly on gaining the optimal control for the unstable power system models and stabilizing them through the Riccati-based feedback stabilization process with sparsity-preserving techniques. We are to find the solution of the Continuous-time Algebraic Riccati Equations (CAREs) governed from the unstable power system models derived from the Brazilian Inter-Connected Power System (BIPS) models, which are large-scale sparse index-1 descriptor systems. We propose the projection-based Rational Krylov Subspace Method (RKSM) for the iterative computation of the solution of the CAREs. The novelties of RKSM are sparsity-preserving computations and the implementation of time-convenient adaptive shift parameters. We modify the Low-Rank Cholesky-Factor integrated Alternating Direction Implicit (LRCF-ADI) technique based nested iterative Kleinman-Newton (KN) method to a sparse form and adjust this to solve the desired CAREs. We compare the results achieved by the Kleinman-Newton method with that of using the RKSM. The applicability and adaptability of the proposed techniques are justified numerically with MATLAB simulations. Transient behaviors of the target models are investigated for comparative analysis through the tabular and graphical approaches.

\end{abstract}

\begin{center}
\keywords{Power system model, Optimal control, Riccati equation, Feedback stabilization, Sparsity preservation}
\end{center}

\section{Introduction}

For the practical purposes, optimal control is a vital part of the engineering interest e.g., industrial control systems, system defense strategy, voltage stability of the power systems, and signal processing \cite{huang2002new,angle2019identifying}. Multi-tasking systems having various components arise in many fields of engineering applications, such as microelectronics, micro-electro-mechanical systems, cybersecurity, computer control of industrial processes, communication systems, etc. These systems are composed of branches of sub-systems and are functioned by very large mathematical models utilizing the interrelated inner mathematical system of higher dimensions.

Power system models are one of the prime branches of the application of optimal controls. For the multi-connected power systems, avoiding cyber-attacks, ensuring the stability of the power connections, and maintaining the compatible frequency level are essential. For those inevitable issues, optimal controls have exigent roles \cite{cui2020deeply}.

The dynamic of a large-scale power system model can be described by the Differential Algebraic Equations (DAEs) as

\begin{equation} \label{eqn:powersystem}
\begin{aligned}

\dot{x}(t) & = f(x_1,x_2,P), \quad f:\mathbb{R}^{n_1 + n_2 + n_3}\rightarrow \mathbb{R}^{n_1},\\
0 & = g(x_1,x_2,P), \quad g:\mathbb{R}^{n_1 + n_2 + n_3}\rightarrow \mathbb{R}^{n_2},

\end{aligned}
\end{equation}

where $x_1\in X_1 \subset \mathbb{R}^{n_1}$ is the vector with differential variables, $x_2\in X_2 \subset \mathbb{R}^{n_2}$ is the vector with algebraic variables and $P \in \mathbb{R}^{n_3}$ is the vector of parameters with $n_1 + n_2 =n$ \cite{kose2003introduction,gross2018steady}. In the state-space representation, the dynamic state variables $x_1$ and instantaneous variables $x_2$ are defined for the specific system, where the parameter $P$ defines the configuration and the operation condition.

The state variables in $x_1$ are time-dependent generator voltages and the parameter $P$ is composed of the system parameters. The control devices together form $ f(x_1,x_2,P)$ and the power flow balance form $g(x_1,x_2,P)$. In case of voltage stability, some equations of the $f(x_1,x_2,P)$ will not be considered. For a fixed parameter $P_0$, linearizing the system (\ref{eqn:powersystem}) around the equilibrium point will provide the following Linear Time-Invariant (LTI) continuous-time system with the input-output equations in the sparse form with the block matrices as

\begin{equation}  \label{eqn:DS_matrix_vector}
\begin{aligned}

\underbrace{ \begin{bmatrix}  E_1  &  0 \\ 0  &  0   \end{bmatrix}}_E

\underbrace{ \begin{bmatrix}  \dot{x_1}(t)  \\ \dot{x_2}(t)    \end{bmatrix}}_{\dot{x}(t)}

&= \underbrace{\begin{bmatrix}  J_1 &  J_2 \\ J_3 & J_4 \end{bmatrix}}_A

\underbrace{ \begin{bmatrix}    x_1(t) \\ x_2(t)  \end{bmatrix}}_{x(t)}
+ \underbrace{\begin{bmatrix}    B_1 \\ B_2  \end{bmatrix}}_B u(t); x(t_0)=x_0,\\

y(t)&= \underbrace{\begin{bmatrix}  C_1  & C_2 \end{bmatrix}}_C

\begin{bmatrix}    x_1(t) \\ x_2(t)  \end{bmatrix} + D u(t),

\end{aligned}

\end{equation}

where $E,A\in \mathbb{R}^{n\times n}$, $B\in \mathbb{R}^{n\times p}$, $C\in \mathbb{R}^{m\times n}$ and $D\in \mathbb{R}^{m\times p}$ with very large $n$ and $p,m\ll n$, represent differential coefficient matrix, state matrix, control multiplier matrix, state multiplier matrix and direct transmission map respectively \cite{mahmud2012full}. In the system (\ref{eqn:DS_matrix_vector}), $x(t)\in \mathbb{R}^n$ is the state vector and $u(t)\in \mathbb{R}^p$ is control (input), while $y(t)\in \mathbb{R}^m$ is the output vector and considering $x(t_0)=x_0$ as the initial state. In most of the state-space representations the direct transmission remains absent and because of that $D=0$. Since $E$ is singular (i.e. $\det(E)=0$), the system (\ref{eqn:DS_matrix_vector}) is called the descriptor system \cite{hossain2019iterative,benner2016reduced}.

Here $x_1\in \mathbb{R}^{n_1}$, $x_2\in \mathbb{R}^{n_2}$ with $n_1+n_2 = n$ are state vectors and other sub-matrices are sparse in appropriate dimensions. If $E_1$ and $J_1$ have full rank, and $J_4$ is non-singular (i.e. $\det(A)\neq 0$), the system is called the index-1 descriptor system \cite{uddin2015computational}. In the current work, we will focus on the stabilization of index-1 descriptor system only. By proper substitution and elimination the descriptor system (\ref{eqn:DS_matrix_vector})  can be converted to the generalized LTI continuous-time system

\begin{equation} \label{eqn:d-state-space}
\begin{aligned}

\mathcal{E}\dot{x}(t)&=\mathcal{A}x(t)+\mathcal{B}u(t),\\
y(t)&=\mathcal{C}x(t)+\mathcal{D}u(t),

\end{aligned}
\end{equation}

where, we have considered the following relations

\begin{equation} \label{eqn:schur_complements}
\begin{aligned}

x:&=x_1, \quad \mathcal{E}:= E_1,\\
\mathcal{A}:&= J_1-J_2{J_4}^{-1}J_3, \quad
\mathcal{B}:= B_1-J_2{J_4}^{-1}B_2,\\
\mathcal{C}:&= C_1-C_2{J_4}^{-1}J_3, \quad
\mathcal{D}:= D-C_2{J_4}^{-1}B_2.

\end{aligned}
\end{equation}

\begin{lemma} [\textbf{Equivalence of transfer functions \cite{gugercin2013model}}] \label{transfer_function}
Assume the transfer functions $G(s)=C(s E - A)^{-1} B + D$  and  $\mathcal{G}(s)=\mathcal{C}(s \mathcal{E} - \mathcal{A})^{-1} \mathcal{B} + \mathcal{D}$ are obtained from the index-1 descriptor system (\ref{eqn:DS_matrix_vector}) and the converted generalized system (\ref{eqn:d-state-space}), respectively. Then, the transfer functions $G(s)$ and $\mathcal{G}(s)$ are identical and hence those systems are equivalent.
\end{lemma}

Though the systems (\ref{eqn:DS_matrix_vector}) and (\ref{eqn:d-state-space}) are equivalent, the system (\ref{eqn:DS_matrix_vector}) is sparse and the system (\ref{eqn:d-state-space}) is dense. This explicit conversion is contradictory with the aim of the work and it will be bypassed by a more efficient way.

LTI continuous-time systems are the pivot ingredient of the present control theory and many other areas of science and engineering \cite{benner2016structure}. Continuous-time Algebraic Riccati Equation (CARE) appears in many branches of engineering applications; especially in electrical fields \cite{chu2011solving,chen2016linear}. The CARE connected to the system (\ref{eqn:d-state-space}) is defined as

\begin{equation} \label{eqn:D-GCARE}
\begin{aligned}

\mathcal{A}^T X\mathcal{E}+\mathcal{E}^T X\mathcal{A}-\mathcal{E}^T X\mathcal{B}\mathcal{B}^T X\mathcal{E}+\mathcal{C}^T \mathcal{C}=0.

\end{aligned}
\end{equation}

If the Hamiltonian matrix corresponding to the system (\ref{eqn:d-state-space}) has no pure imaginary eigenvalues, then the solution $X$ of the CARE (\ref{eqn:D-GCARE}) exists and unique \cite{abou2012matrix}. The solution $X$ of (\ref{eqn:D-GCARE}) is symmetric positive-definite and called stabilizing for the stable closed-loop matrix $\mathcal{A}- (\mathcal{B}\mathcal{B}^T)X\mathcal{E}$. Riccati-based feedback matrix has a prime role in the stabilization approaches for unstable systems \cite{uddin2019riccati,uddin2019efficient}. To find an optimal feedback matrix $K^o$, the Linear Quadratic Regulator (LQR) problem technique can be applied, where the cost functional is defined as

\begin{equation} \label{eqn:LQR}
\begin{aligned}

J(u,x_0)=\int_{0}^{\infty } (x^T(t) \mathcal{C}^T \mathcal{C} x(t)+u^T(t)u(t))dt.

\end{aligned}
\end{equation}

The cost functional (\ref{eqn:LQR}) can be optimized as $J(u^o,x_0)= {x_0}^T Xx_0$ by applying the optimal control $u^o = -K^o x(t)$ generated by the optimal feedback matrix $K^o = \mathcal{B}^T X\mathcal{E}$ associated with the solution matrix $X$ of the CARE (\ref{eqn:D-GCARE}). Using the optimal feedback matrix $K^o$, an unstable LTI continuous-time system can be optimally stabilized by replacing $\mathcal{A}$ by $\mathcal{A}_s = \mathcal{A}-\mathcal{B}K^o$. The stabilized system can be written as

\begin{equation} \label{eqn:st-state-space}
\begin{aligned}

\mathcal{E}\dot{x}(t) = \mathcal{A}_s x(t)+\mathcal{B}u(t), \\
y(t) = \mathcal{C}x(t)+\mathcal{D}u(t).

\end{aligned}
\end{equation}

The eligibility of Rational Krylov Subspace Method (RKSM) for the large-scale LTI continuous-time systems has discussed by Simoncini in \cite{simoncini2007new} and the application of adaptive RKSM to solve large-scale CARE for finding optimal control of the LTI systems has narrated by Druskin \textit{et al.} in \cite{druskin2011adaptive}. Analysis of the basic properties of RKSM for solving large-scale CAREs subject to LTI systems investigated by Simoncini in \cite{simoncini2016analysis}, where the author briefed a new concept of shift parameters that is efficient for the perturbed systems. Very detailed discussion on the numerical solution of large-scale CAREs and LQR based optimal control problems are given by Benner \textit{et al.} in \cite{benner2008numerical}, where the authors narrated the extensions Newton method by means of Alternating Direction Implicit (ADI) technique with convergence properties and the comparative analysis with the RKSM approach. For solving large-scale matrix equations, the Kleinman-Newton method based on Low-Rank Cholesky-Factor ADI (LRCF-ADI) technique is discussed by Kuerschner in \cite{kurschner2016efficient}.

Since there is a scarcity of efficient computational solvers or feasible simulation tools for large-scale CAREs governed from the unstable power system models, the concentration of this work is to develop sparsity-preserving efficient techniques to find the solution of those CAREs. We are proposing a sparsity-preserving and rapid convergent form of the RKSM algorithm for finding solutions of CAREs associated with the unstable power system models and implementation of Riccati-based feedback stabilization. Also, a modified sparse form of LRCF-ADI based Kleinman-Newton method for solving CAREs subject to unstable power system models and corresponding stabilization approach is proposed. The proposed techniques are applied for the stabilization of the transient behaviors of unstable Brazilian Inter-Connected Power System (BIPS) models \cite{varricchio2015hybrid}. Moreover, the comparison of the computational results will be provided in both tabular and graphical methods.

\section{Preliminaries}\label{pre}

In this section, we discuss the background of the proposed techniques, which are derived for generalized LTI continuous-time systems. It includes the derivation of the Rational Krylov Subspace Method (RKSM) technique for solving the Riccati equation, basic structure of the Kleinman-Newton method, and derivation of Low-Rank Cholesky-Factor Alternative Direction Implicit (LRCF-ADI) method for solving Lyapunov equations.

\subsection{Generalized RKSM technique for solving Riccati equations}

In \cite{simoncini2016analysis} Simoncini applied RKSM approach for solving the CARE in the form

\begin{equation} \label{eqn:GCARE}
\begin{aligned}

A^T XE + E^T XA - E^T XBB^T XE + C^T C = 0,

\end{aligned}
\end{equation}

associated with the LTI continuous-time system

\begin{equation} \label{eqn:g-state-space}
\begin{aligned}

E\dot{x}(t) = Ax(t)+Bu(t), \\
y(t) = Cx(t)+Du(t).

\end{aligned}
\end{equation}

If the eigenvalues of the matrix pair $(A,E)$ satisfy $\lambda_i +\bar{\lambda}_j \neq 0, \forall  i,j=1,2, \cdots ,m$, that ensures the solution $X$ of the CARE (\ref{eqn:GCARE}) exists and unique. The orthogonal projector $V\in \mathbb{R}^{n \times m}$ spanned by the $m$-dimensional rational Krylov subspace for a set of shift parameters $\mu_i\in \mathbb{C}^+; i=1,2, \cdots ,m$ is defined as

\small{\begin{align*}

\mathcal{K}_m=span \left(C^T, {(A^T-\mu_1 E^T)}^{-1}C^T,\cdots, \prod_{i=1}^{m} {(A^T-\mu_i E^T)}^{-1} C^T \right).

\end{align*} }

If  $\theta_j$ are the eigenvalues of $(V^T AV, V^T EV)$; $\mathbb{S}_m \in \mathbb{C}^+$ approximates the mirror eigenspace of $A-B(B^T XE)$ and $\delta \mathbb{S}_m$ is its border, the shifts are computed from

\begin{align*}

\mu_{i+1} = arg \left({\mathop{\mathrm{max}}_{\mu \in \delta \mathbb{S}_m} \left|\frac{\prod_{j=1}^i {(\mu-\mu_j)}}{\prod_{j=1}^i {(\mu-\theta_j)}} \right|}\right).

\end{align*}

According to the Galerkin condition and after simplification by matrix algebra, a low-rank CARE can be obtained as

\begin{equation} \label{eqn:LR-GCARE}
\begin{aligned}

\hat{A}^T \hat{X}\hat{E} +\hat{E}^T \hat{X} \hat{A}-\hat{E}^T \hat{X} \hat{B}\hat{B}^T \hat{X}\hat{E}+ \hat{C}^T \hat{C}=0,

\end{aligned}
\end{equation}

where $\hat{X}=V^T XV, \hat{E}=V^T EV, \hat{A}=V^T AV, \hat{B}=V^T B$ and $\hat{C}=CV$. The equation (\ref{eqn:LR-GCARE}) is an approximated low-rank CARE and can be solved by any conventional method or MATLAB \texttt{care} command. Here $\hat{X}\in \mathbb{R}^{m \times m}$ is taken as low-rank approximation of $X$, corresponding to the low-rank CARE (\ref{eqn:LR-GCARE}). Then residual of the $m$-th iteration is

\begin{align*}

\|{\mathcal{R}_m}\|_F = \|{SJS^T}\|_F; \quad J&= \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix},

\end{align*}

where ${\|{.}\|}_F$ denotes the Frobenius norm and $S$ is a block upper triangular matrix in the QR factorization of the matrix $U$ derived as

\begin{align*}

U ={ \begin{bmatrix} v_{m+1} \mu_{m+1} \\  E^T V_m \hat{X} H_m^{-T} e_m h_{m+1,m}^T  \\  -(I_n - V_m V_m^T) A^T v_{m+1}
\end{bmatrix} }^T,

\end{align*}

where $H_m$ is a block upper Hessenberg matrix and $e_m$ is the matrix formed by the last $p$ columns of the $mp$-order identity matrix. For  $C^T = Q_0 R_0$ such that $R_0 = \beta_0$, the relative-residual can be estimated as

\begin{align*}

{\|{\mathcal{R}_m}\|}_{(relative)} = \frac{\|{\mathcal{R}_m\|}_F}{\|{\beta_0^T \beta_0}\|_F}.

\end{align*}

Through RKSM the low-rank factor $Z$ of the approximate solution $X$ of the CARE (\ref{eqn:GCARE}) needs to be estimated, such that $X \approx Z Z^T $. The low-rank solution $\hat{X}$ is symmetric positive definite and the original solution $X$ can be approximated as $X=V \hat{X} V^T$. By the eigenvalue decomposition to the approximate solution $\hat{X}$ and truncating the negligible eigenvalues, the possible lowest order factor $Z$ of $X$ can be estimated as

\begin{align*}

X &= V \hat{X} V^T = V (T\Lambda T^T) V^T \\ &= V \begin{bmatrix} T_1 & T_2 \end{bmatrix} \begin{bmatrix} \Lambda_1 & 0 \\ 0 & \Lambda_2\end{bmatrix} \begin{bmatrix} T_1^T \\ T_2^T \end{bmatrix} V^T \\ &= VT_1 \Lambda_1 T_1^T V^T =  (VT_1 \Lambda_1^{\frac{1}{2}})(VT_1 \Lambda_1^{\frac{1}{2}})^T \\ &= ZZ^T.

\end{align*}

Here, $\Lambda_2$ consists the negligible eigenvalues. Summary of the above process is given in the Algorithm-\ref{pre:RKSM.}.

\begin{algorithm}[t]
\SetAlgoLined
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\caption{Generalized RKSM.}
\label{pre:RKSM.}
\Input{ $E, A, B, C, i_{max}$ (number of iterations) and $\mu_i$ (initial shifts).}
\Output{Low-rank factored solution $Z$ such that $X\approx ZZ^T$.}

Compute $Q_0 R_0 = C^T$ (QR factorization).\\
Choose $V_0 = Q_0$.\\

\While{{not converged \ or \ $m \leq i_{\text{max}}$}}{
Solve $v_{m+1} = {(A^T-\mu_{m+1} E^T)}^{-1}V_m$.\\
Compute shift for the next iteration.\\
Using Arnoldi algorithm orthogonalize $v_{m+1}$ against $V_m$ to obtain $\hat{v}_{m+1}$, such that $V_{m+1}=\begin{bmatrix} V_m, \hat{v}_{m+1} \end{bmatrix}$.\\
Assuming $\hat{E}=V_{m+1}^T EV_{m+1}$, $\hat{A}=V_{m+1}^T AV_{m+1}$, $\hat{B}=V_{m+1}^TB$ and $\hat{C}=CV_{m+1}$, for $\hat{X}$ solve the reduced-order Riccati equation $\hat{A}^T \hat{X} \hat{E}+ \hat{E}^T \hat{X} \hat{A}- \hat{E}^T \hat{X} \hat{B} \hat{B}^T \hat{X} \hat{E} + \hat{C}^T \hat{C}=0$.\\
Compute ${\|{\mathcal{R}_m}\|}_{(\text{relative})}$ for convergence.\\
}

Compute eigenvalue decomposition $\hat{X} = T\Lambda T^T=\begin{bmatrix} T_1 & T_2 \end{bmatrix} \begin{bmatrix} \Lambda_1 & 0 \\ 0 & \Lambda_2\end{bmatrix} \begin{bmatrix} T_1^T \\ T_2^T \end{bmatrix}$.\\
For negligible eigenvalues truncate $\Lambda_2$ and compute $Z=V_{m+1} T_1 \Lambda_1^{\frac{1}{2}}$. \\
\end{algorithm}

\subsection{Modified LRCF-ADI method for solving Lyapunov equations}

The generalized Continuous-time Algebraic Lyapunov Equation (CALE) associated with the LTI continuous-time system (\ref{eqn:g-state-space}) is

\begin{equation} \label{eqn:GCALE}
\begin{aligned}

A^T XE+ E^TXA =- C^T C.

\end{aligned}
\end{equation}

The shift parameters $\mu_i \in \mathbb{C}^-$ are allowed and the initial iteration is taken as $X_0={X_0}^T \in \mathbb{R}^{n\times n}$. Assume $Z_i \in \mathbb{R}^{n\times ip}$ as the low-rank Cholesky factor of $X_i \in \mathbb{R}^{n\times n}$ such that $X_i=Z_i Z_i^T$ \cite{li2000model}. The ADI algorithm can be formed in terms of Cholesky factor $Z_i$ of $X_i$ and there will be no need to estimate or store $X_i$ at each iteration as only $Z_i$ is required \cite{penzl1999cyclic}.

Considering $\gamma_i =\sqrt{-2\Rel(\mu_i)}$, the conventional low-rank Cholesky factor ADI iterations yield the form as

\begin{equation} \label{eqn:ADI-factor4}
\begin{aligned}

V_1 &= (A^T+\mu_1 E^T)^{-1}C^T,\\
Z_1 &=\gamma_1 V_1=\sqrt{-2\Rel(\mu_i)} (A^T+\mu_1 E^T)^{-1}C^T,\\
V_i &=H_{i-1,i} V_{i-1}=V_{i-1}-(\mu_i + \bar{\mu}_{i-1})(A^T+\mu_i E^T)^{-1} E^T V_{i-1},\\
Z_i &=\begin{bmatrix} Z_{i-1} \quad \gamma_i V_i \end{bmatrix}=\begin{bmatrix} Z_{i-1} \quad \sqrt{-2\Rel(\mu_i)} V_i \end{bmatrix}.

\end{aligned}
\end{equation}

Assume a set of adjustable shift parameters, for two subsequent block iterates $V_i, V_{i+1}$  of the ADI technique related to the pair of complex conjugated shifts $\{ \mu_i, \mu_{i+1}:=\bar{\mu}_i\}$ it holds

\begin{equation} \label{eqn:t1}
\begin{aligned}

V_{i+1}=\bar{V}_i + 2\delta_i \text{Im}(V_i),

\end{aligned}
\end{equation}

where $\bar{(.)}$ indicates the complex conjugate with $\delta_i=\frac{\text{Re}(\mu_i)}{\text{Im}(\mu_i)}$ and iterates associated to real shifts are always purely real \cite{benner2013efficient}.

Then the following matrix for the basis extension can be obtained

\begin{equation} \label{eqn:iterations}
\begin{aligned}

\begin{bmatrix} V_i \ V_{i+1} \end{bmatrix}=\begin{bmatrix} \gamma_i (\Rel(V_i)+\delta_i \Imag(V_i)) \quad \gamma_i \sqrt{\delta_i^2 + 1} \Imag(V_i)\end{bmatrix}.

\end{aligned}
\end{equation}

Then, for a pair of complex conjugate shifts at any iteration, the low-rank factor $Z_i$ can be computed as

\begin{equation} \label{eqn:ADI-factor5}
\begin{aligned}

Z_{i+1}=\begin{bmatrix}Z_{i-1} \ \gamma_i (\Rel(V_i)+\delta_i \Imag(V_i)) \quad \gamma_i \sqrt{\delta_i^2 + 1} \Imag(V_i)\end{bmatrix}.

\end{aligned}
\end{equation}

The modified techniques discussed above is summarized in Algorithm-(\ref{pre:G-LRCF-ADI2}).

\begin{algorithm}[t]
\SetAlgoLined
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\caption{Modified LRCF-ADI.}
\label{pre:G-LRCF-ADI2}
\Input{ $E, A, C,  i_{max}$ (number of iterations) and shift parameters $\{\mu_j\}_{j=1}^{i_{max}}$.}
\Output{Low-rank Cholesky-factor $Z$ such that $X\approx ZZ^T$.}

Consider $Z_0=[ \ ]$.\\

\For {$i \gets 1$ to $i_{max}$} {
\eIf{$i=1$}{Solve $V_1=(A^T+\mu_1 E^T)^{-1} C^T$.}
{Compute $V_i=V_{i-1}-(\mu_i + \bar{\mu}_{i-1})(A^T+\mu_i E^T)^{-1} E^T V_{i-1}$.}

\eIf{$\text{Im}(\mu_i)=0$}{Update $Z_i=\begin{bmatrix} Z_{i-1} \quad \sqrt{-2\mu_i} V_i \end{bmatrix}$.}

{Assume $\gamma_i =\sqrt{-2\text{Re}(\mu_i)}$,  $\delta_i=\frac{\text{Re}(\mu_i)}{\text{Im}(\mu_i)}$,\\
Update $Z_{i+1}=\begin{bmatrix}Z_{i-1} \quad \gamma_i (\text{Re}(V_i)+\delta_i \text{Im}(V_i)) \quad \gamma_i \sqrt{\delta_i^2 + 1} \text{Im}(V_i) \end{bmatrix}$,\\
Compute $V_{i+1}=\bar{V}_i + 2\delta_i \ \text{Im}(V_i)$.\\
$i=i+1$
}

}

\end{algorithm}

\subsection{Basic Kleinman-Newton method}

For the system (\ref{eqn:g-state-space}), the residual $\mathcal{R}(X)$ of the CARE (\ref{eqn:GCARE}) is

\begin{equation} \label{eqn:n-res}
\begin{aligned}

\mathcal{R}(X)=A^T X E + E^T XA - E^T XBB^T X E + C^T C.

\end{aligned}
\end{equation}

Form the Fre'chet derivative at $X=X_i$, we have

\begin{equation} \label{eqn:def-res}
\begin{aligned}

\mathcal{R}'(X_i)=(A-BB^T X E)^T X_i E + E^T X_i (A-BB^T X E).

\end{aligned}
\end{equation}

Consider the Newton iteration and apply (\ref{eqn:def-res}), we have

\begin{equation} \label{eqn:newton}
\begin{aligned}

\mathcal{R}'(\Delta X_i)+\mathcal{R}(X_i)=0.

\end{aligned}
\end{equation}

Now, put $\Delta X_i = X_{i+1} - X_i$ in (\ref{eqn:newton}). Then, after simplification we get

\begin{equation} \label{eqn:n-care}
\begin{aligned}

& (A-BB^T X_i E)^T X_{i+1} E + E^T X_{i+1} (A-BB^T X_i E) \\ &= -C^T C-E^T X_i BB^T X_i E.

\end{aligned}
\end{equation}

Then, assume $\tilde{A}_i=A-BB^T X_i E$ and $W_i= \begin{bmatrix}  C^T \quad  E^T X_i B \end{bmatrix}$, then equation (\ref{eqn:n-care}) reduces to a generalized Lyapunov equation such as

\begin{equation} \label{eqn:n-lyap}
\begin{aligned}

\tilde{A}_i^T X_{i+1} E + E^T X_{i+1} \tilde{A}_i = - W_i W_i^T.

\end{aligned}
\end{equation}

The generalized CALE (\ref{eqn:n-lyap}) can be solved for $X_{i+1}$ by any conventional method, such as Low-Rank Cholesky-Factor Alternative Direction Implicit (LRCF-ADI) method and the corresponding feedback matrix $K_{i+1}=B^T X_{i+1}E$ can be estimated. The whole mechanism is called the Kleinman-Newton method \cite{benner2014numerical} for solving generalized CARE. The summary of the method is given in Algorithm-\ref{pre:KN.}.

\begin{algorithm}[t]
\SetAlgoLined
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\caption{Basic Kleinman-Newton.}
\label{pre:KN.}
\Input{$E,A,B,C$ and $X_0$ (initial assumption).}
\Output{Approximate solution $X$ and feedback matrix $K$.}

\While{ $i \leq i_{\text{max}} $}{
Compute $\tilde{A}_i=A-BB^T X_i E$ and $W_i= \begin{bmatrix} C^T \quad     E^T X_i B \end{bmatrix}$;\\
For $X_{i+1}$, solve $\tilde{A}_i^T X_{i+1} E + E^T X_{i+1} \tilde{A}_i = - W_i W_i^T$;\\
Compute $K_{i+1}=B^T X_{i+1} E$.
}
\end{algorithm}

\section{Solving Riccati equations arising from the index-1 descriptor systems}

In this section, we discuss the updated RKSM techniques for solving the Riccati equation derived from index-1 descriptor systems, it includes the stopping criteria, sparsity preservation, and estimation of the optimal feedback matrix. Then, LRCF-ADI based Kleinman-Newton method with the adjustment for index-1 descriptor systems and finding the optimal feedback matrix are discussed. Finally, the stabilized system is narrated accordingly.

\subsection{Updated Rational Krylov subspace method}

Let us consider the LTI continuous-time system (\ref{eqn:d-state-space}) and introduce an orthogonal projector $V$ spanned by the $m$ dimensional rational Krylov subspace for a set of given shift parameters $\mu_i\in \mathbb{C}^+; i=1,2,...,m$ is defined as

\small{\begin{align*}

\mathcal{K}_m=span \left(\mathcal{C}^T, {(\mathcal{A}^T-\mu_1 \mathcal{E}^T)}^{-1}\mathcal{C}^T,...,\prod_{i=1}^{m} {(\mathcal{A}^T-\mu_i \mathcal{E}^T)}^{-1} \mathcal{C}^T \right).

\end{align*} }

Again, consider the CARE (\ref{eqn:D-GCARE}) and apply the Galerkin condition on it. Then, after the simplification by matrix algebra, a low-rank CARE can be achieved as

\begin{equation} \label{eqn:LR-DCARE}
\begin{aligned}

\hat{\mathcal{A}}^T \hat{X} \hat{\mathcal{E}}+\hat{\mathcal{E}}^T \hat{X} \hat{\mathcal{A}}-\hat{\mathcal{E}}^T \hat{X} \hat{\mathcal{B}} \hat{\mathcal{B}}^T \hat{X} \hat{\mathcal{E}}+\hat{\mathcal{C}}^T \hat{\mathcal{C}}=0,

\end{aligned}
\end{equation}

where, $\hat{X}=V^T XV, \hat{\mathcal{E}}=V^T \mathcal{E}V, \hat{\mathcal{A}}=V^T \mathcal{A}V, \hat{\mathcal{B}}=V^T \mathcal{B}$ and $\hat{\mathcal{C}}=\mathcal{C}V$. The equation (\ref{eqn:LR-DCARE}) is a low-rank CARE and can be solved by MATLAB \texttt{care} command or any existing methods, such as \texttt{Schur}-decomposition method.

For the quick and smooth convergence of the proposed algorithm, adjustable shift selection is crucial and we are adopting the adaptive shift approach for index-1 descriptor systems \cite{benner2014self}. This process required to be recursive and in each step the subspace to all the projector $V$ generated with the current set of shifts will be extended.

\subsubsection{Stopping criteria and related theorem}

Arnoldi relation is a very essential tool for the computation of residual of the RKSM iterations. To avoid extra matrix-vector multiplies per iteration, the computation of projected matrix $T_m$ can be performed more efficiently than the explicit product $T_m = V_m^T \mathcal{A}V_m$. To find the stopping criteria, we have to consider following lemmas.

\begin{lemma} [\textbf{Arnoldi relation \cite{abidi2016adaptive}}] \label{a.relation}
Let $\mathcal{K}_m = \text{span} (V_m)$ be the rational Krylov subspace with the shift parameters $\mu_i\in \mathbb{C}^+; i=1,2, \cdots ,m$. Then $\mathcal{K}_m$ satisfies the Arnoldi relation as follows

\begin{equation} \label{eqn:a.relation1}
\begin{aligned}

\mathcal{A}^TV_m=V_m T_m + \hat{v}_{m+1} g_m^T \ ; \quad V_m^T V_m = I_m,

\end{aligned}
\end{equation}

where $\hat{v}_{m+1} \beta_1 = v_{m+1} \mu_{m+1} - (I_n - V_m V_m^T) \mathcal{A}^T v_{m+1}$ is the QR decomposition of the right hand side matrix with $g_m^T = \beta_1 h_{m+1,m} e_m^T H_m^{-1}$.
\end{lemma}

\begin{lemma} [\textbf{Building the projected matrix \cite{lin2013minimal}}] \label{projected}
Let the column vectors of $V_m$ be an orthonormal basis of the rational Krylov subspace $\mathcal{K}_m$ with the block diagonal matrix $D_m = \text{diag} (\mu_1 \mathcal{E}^T, \mu_2 \mathcal{E}^T, \cdots , \mu_m \mathcal{E}^T)$, where $\{ \mu_1, \mu_2, \cdots , \mu_m \}$ is the set of shift parameters used in the algorithm. Then for the projected matrix $T_m = V_m^T \mathcal{A} V_m$ the following relation holds

\begin{equation} \label{eqn:projected1}
\begin{aligned}

T_m = \left[ I_m + V_m^T \mathcal{E}^T V_m H_m D_m - V_m^T \mathcal{A}^T v_{m+1} h_{m+1,m} e_m^T \right] H_m^{-1}.

\end{aligned}
\end{equation}

\end{lemma}

\begin{theorem}  [\textbf{Residual of the RKSM iterations}] \label{Trksm}

Let $V_m$ be the orthogonal projector spanned by the rational Krylov subspace $\mathcal{K}_m$ and $X\approx V\hat{X} V^T$ is the solution of the CARE (\ref{eqn:D-GCARE}) using the low-rank solution $\hat{X}$. Then, the residual of $m$-th iteration can be computed as

\begin{equation} \label{eqn:trksm1}
\begin{aligned}

{\|{\mathcal{R}_m}\|}_F={\|{ S J S^T}\|}_F \ ; \quad J=\begin{bmatrix} 0 & I & 0 \\ I & 0 & I \\  0 & I & 0 \end{bmatrix},

\end{aligned}
\end{equation}

where ${\|{.}\|}_F$ denotes the Frobenius norm and $S$ is a block upper triangular matrix in the QR factorization of the matrix $U$ is defined as

\begin{equation} \label{eqn:trksm2}
\begin{aligned}

U = \begin{bmatrix} v_{m+1} \mu_{m+1} \quad \mathcal{E}^T V_m \hat{X} H_m^{-T} e_m h_{m+1,m}^T \quad -(I_n - V_m V_m^T) \mathcal{A}^T v_{m+1} \end{bmatrix}.

\end{aligned}
\end{equation}

\end{theorem}

\begin{proof}

Assume $f = (I_n - V_m V_m^T) \mathcal{A}^T v_{m+1}$ and consider the reduced QR factorization $\mathcal{C}^T=V_0 \beta_0$. Then, by putting the relations in equation (\ref{eqn:a.relation1}) of Lemma-(\ref{a.relation}), the relation can be written as

\begin{equation} \label{eqn:trksm3}
\begin{aligned}

\mathcal{A}^TV_m &= V_m T_m + {v}_{m+1} \mu_{m+1} h_{m+1,m} e_m^T H_m^{-1} -  g h_{m+1,m} e_m^T H_m^{-1}\\
&= V_m T_m + ({v}_{m+1} \mu_{m+1} - f) h_{m+1,m} e_m^T H_m^{-1}.

\end{aligned}
\end{equation}

The residual of the CARE (\ref{eqn:D-GCARE}) can be written as

\begin{equation} \label{eqn:trksm4}
\begin{aligned}

\mathcal{R}(X) = \mathcal{A}^T X\mathcal{E} + \mathcal{E}^T X\mathcal{A} - \mathcal{E}^T X\mathcal{B} \mathcal{B}^T X\mathcal{E} + \mathcal{C}^T \mathcal{C}.

\end{aligned}
\end{equation}

Consider the approximate solution using the low-rank solution $\hat{X}$ as $X= V_m \hat{X} V_m^T$ and equation-(\ref{eqn:projected1}) in Lemma-(\ref{projected}), then applying (\ref{eqn:trksm3}) in (\ref{eqn:trksm4}), we get

\begin{equation} \label{eqn:trksm5}
\begin{aligned}

&\mathcal{R}(X_m) = V_m T_m \hat{X}_m V_m^T \mathcal{E} + ({v}_{m+1} \mu_{m+1} - f) h_{m+1,m} e_m^T H_m^{-1} \hat{X}_m V_m^T \mathcal{E} + \mathcal{E}^T V_m \hat{X} T_m^T V_m^T \\ &+ \mathcal{E}^T V_m \hat{X} H_m^{-T} e_m h_{m+1,m}^T ({v}_{m+1} \mu_{m+1} - f)^T - \mathcal{E}^T V_m \hat{X} V_m^T \mathcal{B} \mathcal{B}^T V_m \hat{X} V_m^T \mathcal{E} + \mathcal{C}^T \mathcal{C}\\
&= ({v}_{m+1} \mu_{m+1} - f) h_{m+1,m} e_m^T H_m^{-1} \hat{X}_m V_m^T \mathcal{E} + \mathcal{E}^T V_m \hat{X} H_m^{-T} e_m h_{m+1,m}^T ({v}_{m+1} \mu_{m+1} - f)^T \\ &+ \mathcal{M}^T X \mathcal{E} + \mathcal{E}^T X \mathcal{M} - \mathcal{E}^T X \mathcal{B} \mathcal{B}^T X \mathcal{E} + \mathcal{C}^T \mathcal{C} \ ; \ \mathcal{M} = \mathcal{A}^T V_m V_m^T\\
&= \mathcal{E}^T V_m \hat{X} H_m^{-T} e_m h_{m+1,m}^T \mu_{m+1}^T v_{m+1}^T + ({v}_{m+1} \mu_{m+1} - f) h_{m+1,m} e_m^T H_m^{-1} \hat{X}_m V_m^T \mathcal{E} \\ &- \mathcal{E}^T V_m \hat{X} H_m^{-T} e_m h_{m+1,m}^T f^T + 0 \\
&=\begin{bmatrix} v_{m+1} \mu_{m+1} \ \ \mathcal{E}^T V_m \hat{X}_m H_m^{-T} e_m h_{m+1,m}^T \ \ -f \end{bmatrix}
\begin{bmatrix} 0 & I & 0 \\ I & 0 & I \\ 0 & I & 0 \end{bmatrix}
\begin{bmatrix} \mu_{m+1}^T v_{m+1}^T \\  h_{m+1,m} e_m^T H_m^{-1} \hat{X}_m^T V_m^T \mathcal{E} \\ -f^T \end{bmatrix}\\
&= S J S^T.

\end{aligned}
\end{equation}

Thus, the proof follows from the Frobenius norm of the equation (\ref{eqn:trksm5}).

\end{proof}

Then the relative-residual with respect to $\beta_0$ can be estimated as follows

\begin{equation} \label{eqn:relres}
\begin{aligned}

{\|{\mathcal{R}_m}\|}_{(\text{relative})} = \frac{\|{\mathcal{R}_m}\|_F}{\|{\beta_0^T \beta_0}\|_F}.

\end{aligned}
\end{equation}

\subsubsection{Sparsity preservation}

The matrix $\mathcal{A}$ in (\ref{eqn:d-state-space}) is in dense form, which is contradictory to the aim of the work and the rate of convergence of the converted system is slow enough. So, to bypass these drawbacks at each iteration a shifted linear system needs to be solved for $v_i$ as

\begin{equation} \label{eqn:sparse_shifted}
\begin{aligned}

(\mathcal{A}^T-\mu_i \mathcal{E}^T)v_i &=V_{i-1},\\
or, {\begin{bmatrix}  J_1-{\mu_i}E_1  &  J_2 \\ J_3  &  J_4   \end{bmatrix}}^T \begin{bmatrix}  v_i  \\ \Gamma \end{bmatrix}
&=  \begin{bmatrix}    V_{i-1} \\ 0  \end{bmatrix}.
\end{aligned}

\end{equation}

Here $\Gamma = - J_4^{-1} J_3 v_i$ is the truncated term. The linear system (\ref{eqn:sparse_shifted}) is higher dimensional but sparse and can be solved by the conventional sparse-direct solvers very efficiently \cite{uddin2019computational}. To improve the consistency of the RKSM approach, explicit form of the reduced-order matrices must not be used to construct reduced-order system. The sparsity-preserving reduced-order matrices can be attained by following way

\begin{equation} \label{eqn:sprom}
\begin{aligned}

\hat{\mathcal{E}} &= V^T E_1 V,\quad
\hat{\mathcal{A}} = V^T J_1 V - (V^T J_2){J_4}^{-1}(J_3 V),\\
\hat{\mathcal{B}} &= V^T B_1 - (V^T J_2) {J_4}^{-1} B_2,\quad
\hat{\mathcal{C}} = C_1 V - C_2 {J_4}^{-1} (J_3 V).

\end{aligned}
\end{equation}

\subsubsection{Treatment for the unstable systems}

If the system is unstable, a Bernoulli stabilization is required through an initial-feedback matrix $K_0$ to estimate $\mathcal{A}_f = \mathcal{A}-\mathcal{B} K_0$ and the matrix $\mathcal{A}$ needs to be replaced \cite{benner2017model}. Then, the system (\ref{eqn:d-state-space}) and CARE (\ref{eqn:D-GCARE}) need to be re-defined as

\begin{equation} \label{eqn:re-state-space}
\begin{aligned}

\mathcal{E}\dot{x}(t)&=\mathcal{A}_f x(t)+\mathcal{B}u(t),\\
y(t)&=\mathcal{C}x(t)+\mathcal{D}u(t),

\end{aligned}
\end{equation}

\begin{equation} \label{eqn:re-GCARE}
\begin{aligned}

\mathcal{A}_f^T X\mathcal{E}+\mathcal{E}^T X\mathcal{A}_f-\mathcal{E}^T X\mathcal{B} \mathcal{B}^T X\mathcal{E}+\mathcal{C}^T \mathcal{C}=0.

\end{aligned}
\end{equation}

Then, for every iterations $K=\hat{\mathcal{B}}^T \hat{X}V^T \mathcal{E}$ needs to be updated by the solution $\hat{X}$ of (\ref{eqn:LR-GCARE}) and the rational Krylov subspace for the projector $V$ needs to be redefined as

\small{\begin{align*}

\mathcal{K}_m=span \left(\mathcal{C}^T, {(\mathcal{A}_f^T-\mu_1 \mathcal{E}^T)}^{-1}\mathcal{C}^T,...,\prod_{i=1}^{m} {(\mathcal{A}_f^T-\mu_i \mathcal{E}^T)}^{-1} \mathcal{C}^T \right).

\end{align*} }

For the stabilized system using the initial-feedback matrix $K_0$, the expressions (\ref{eqn:sparse_shifted}) can be written as

\begin{equation} \label{eqn:sparse_shifted_st}
\begin{aligned}

(\mathcal{A}_f^T-\mu_i \mathcal{E}^T)v_i &=V_{i-1},\\
or, {\begin{bmatrix}  (J_1-B_1 K_0)-{\mu_i}E_1  &  J_2 \\ J_3-B_2 K_0  &  J_4   \end{bmatrix}}^T \begin{bmatrix}  v_i  \\ *\end{bmatrix} &=  \begin{bmatrix} V_{i-1} \\ 0  \end{bmatrix}.

\end{aligned}
\end{equation}

To evaluate the shifted linear systems, explicit inversion of $\mathcal{A} - \mathcal{B} K$ should be avoided in practice, instead the \textit{Sherman-Morrison-Woodbury} formula needs to be used as follows

\begin{align*}
(\mathcal{A} - \mathcal{B} K)^{-1} = \mathcal{A}^{-1} + \mathcal{A}^{-1} \mathcal{B}(I- K\mathcal{A}^{-1} \mathcal{B})^{-1}K \mathcal{A}^{-1}.
\end{align*}

\subsubsection{Estimation of the optimal feedback matrix}

The low-rank solution $\hat{X}$ is symmetric and positive definite and can be factorized  as $\hat{X}=YY^T$. The original solution can be approximated as $X=V \hat{X} V^T=VY(VY)^T$. Finally, the desired low-rank factored solution $Z=VY$ of the CARE (\ref{eqn:D-GCARE}) will be stored and the optimal feedback matrix $K^o = \mathcal{B}^T X \mathcal{E} = \mathcal{B}^T (ZZ^T) \mathcal{E}$ can be estimated. This process is iterative and will continue until the desired convergence is achieved. The whole process is summarized in the Algorithm-\ref{pre:RKSM-FB.}.

\begin{algorithm}[t]
\SetAlgoLined
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\caption{Updated RKSM (sparsity-preserving).}
\label{pre:RKSM-FB.}
\Input{ $E_1, J_1, J_2, J_3, J_4, B_1, B_2, C_1, C_2, K_0$ (initial feedback matrix) $i_{max}$ (number of iterations) and $\mu_i$ (initial shifts).}
\Output{Low-rank factored solution $Z$ such that $X\approx ZZ^T$ and optimal feedback matrix $K^o$.}

Compute $Q_0 R_0 = (C_1 - C_2{J_4}^{-1}J_3)^T$ (QR factorization).\\
Choose $V_0 = Q_0$.\\
Choose $K=K_0$.\\

\While{{not converged \ or \ $m \leq i_{\text{max}}$}}{
Solve the linear system (\ref{eqn:sparse_shifted_st}) for $v_{m+1}$.\\
Compute adaptive shifts for the next iterations (if store is empty).\\
Using Arnoldi algorithm orthogonalize $v_{m+1}$ against $V_m$ to obtain $\hat{v}_{m+1}$, such that $V_{m+1}=\begin{bmatrix} V_m, \hat{v}_{m+1} \end{bmatrix}$.\\
Assuming $\hat{\mathcal{E}}$, $\hat{\mathcal{A}}$, $\hat{\mathcal{B}}$ and $\hat{\mathcal{C}}$ are defined in (\ref{eqn:sprom}), for $\hat{X}$ solve the reduced-order Riccati equation
$\hat{\mathcal{A}}^T \hat{X}\hat{\mathcal{E}}+ \hat{\mathcal{E}}^T \hat{X} \hat{\mathcal{A}}-\hat{\mathcal{E}}^T \hat{X} \hat{\mathcal{B}} \hat{\mathcal{B}}^T \hat{X} \hat{\mathcal{E}} + \hat{\mathcal{C}}^T \hat{\mathcal{C}}=0$.\\
Update $K=(\hat{\mathcal{B}}^T \hat{X})V_{m+1} E_1$.\\
Compute ${\|{\mathcal{R}_m}\|}_{(\text{relative})}$ for convergence.\\
}
Compute eigenvalue decomposition
$\hat{X} = T\Lambda T^T=\begin{bmatrix} T_1 & T_2 \end{bmatrix}
\begin{bmatrix} \Lambda_1 & 0 \\ 0 & \Lambda_2\end{bmatrix}
\begin{bmatrix} T_1^T \\ T_2^T \end{bmatrix}$.\\
For negligible eigenvalues truncate $\Lambda_2$ and construct
$Z=V_{m+1} T_1 \Lambda_1^{\frac{1}{2}}$.\\
Compute the optimal feedback matrix $K^o = (B_1 - J_2{J_4}^{-1}B_2)^T(ZZ^T)E_1$.
\end{algorithm}

\subsection{Updated LRCF-ADI Based modified Kleinman-Newton method}

In each iteration of Algorithm-\ref{pre:KN.}, the generalized CALE needs to be solved for once and there are several techniques available to do it. In his Ph.D. thesis, Kuerschner discussed low-rank ADI approaches (Algorithm-3.2 chapter-3 and Algorithm-6.2 in chapter-6) for solving generalized CALE derived from generalized CARE in the iterative loops of Kleinman-Newton algorithm \cite{kurschner2016efficient}. Now, we need to implement above mechanisms for the index-1 descriptor system in the sparse form. For the adjustment, some modifications are required as given below.

\subsection{Convergence criteria and recurrence relations}

The computation of residuals of the ADI iterations can be achieved using the simplified technique implementing the adjustable shift parameters. This approach will be efficient for time-dealing and memory allocation. Following lemmas represent some important properties of the ADI iterates.

\begin{lemma} [\textbf{Convergence of ADI iterates \cite{feitzinger2009inexact}}] \label{Lres1}
Let $X_i$ be the solution of the CALE (\ref{eqn:GCALE}) and consider $X_i^{(k)}$ be an iterate of the ADI method. Then for all $\mu_i \in \mathbb{C}^-$ the relation holds

\begin{equation} \label{eqn:lyapres1}
\begin{aligned}

X_i^{(k+1)} - X_i = \left( \prod_{j=1}^{k}{\mathcal{A}_{k,\mu_j}} \right) (X_l^{(0)} - X_l) \left( \prod_{j=1}^{k}{\mathcal{A}_{k,\mu_j}} \right)^T,

\end{aligned}
\end{equation}

where $\mathcal{A}_{k,\mu_j} = (\mathcal{A}^T-\bar{\mu_j} \mathcal{E}^T)(\mathcal{A}^T+\mu_j \mathcal{E}^T)^{-1}$.
\end{lemma}

\begin{lemma} [\textbf{Residual of the ADI iterates \cite{feitzinger2009inexact}}] \label{Lres2}
Let $X_i^{(k)}$ be an iterate of the CALE (\ref{eqn:GCALE}) by the ADI method. Then considering $\mathcal{A}_{k,\mu_j} = (\mathcal{A}^T-\bar{\mu_j} \mathcal{E}^T)(\mathcal{A}^T+\mu_j \mathcal{E}^T)^{-1}$ for all $\mu_i \in \mathbb{C}^-$ the residuals at $X_i^{(k)}$ have the form

\begin{equation} \label{eqn:lyapres2}
\begin{aligned}

\mathcal{R}_i^{(k+1)} &= \mathcal{E}^T X_i^{(k)} \mathcal{A} + \mathcal{A}^T X_i^{(k)} \mathcal{E} + \mathcal{C}^T \mathcal{C} \\
&= \left( \prod_{j=1}^{k}{\mathcal{A}_{k,\mu_j}} \right) (\mathcal{E}^T X_i^{(0)} \mathcal{A} + \mathcal{A}^T X_i^{(0)} \mathcal{E} + \mathcal{C}^T \mathcal{C}) \left( \prod_{j=1}^{k}{\mathcal{A}_{k,\mu_j}} \right)^T.

\end{aligned}
\end{equation}
\end{lemma}

\begin{theorem} [\textbf{Residual factor of the ADI iterations}] \label{Tadi2}
The residual for CALE (\ref{eqn:GCALE}) at $i$-th iteration of ADI method is of rank at most $m$ and is given by

\begin{equation} \label{eqn:residue1}
\begin{aligned}

\mathcal{R}(X_i) = \mathcal{A}^T X_i \mathcal{E}+\mathcal{E}^T X_i \mathcal{A}+\mathcal{C}^T \mathcal{C} = \mathcal{W}_i \mathcal{W}_i^T,

\end{aligned}
\end{equation}

with $\mathcal{W}_i=\left( \prod_{i=1}^{m} {(\mathcal{A}^T-\bar{\mu_i} \mathcal{E}^T)(\mathcal{A}^T+\mu_i \mathcal{E}^T)^{-1}} \right)\mathcal{C}^T$, for all $\mu_i \in \mathbb{C}^-$.
\end{theorem}

\begin{proof}
Consider $X$ is the solution of the CALE (\ref{eqn:GCALE}) and $X_i$ is its $i$-th iterate by ADI method. Then the residual of the $i$-th iteration in terms of $X_i$ can be written as

\begin{equation} \label{eqn:residue2}
\begin{aligned}

\mathcal{R}(X_i) = \mathcal{A}^T X_i \mathcal{E} + \mathcal{E}^T X_i \mathcal{A} + \mathcal{C}^T \mathcal{C}.

\end{aligned}
\end{equation}

For all $\mu_i \in \mathbb{C}^-$, the Stein's equation is equivalent to the CALE (\ref{eqn:GCALE}) and it can be written as

\begin{equation} \label{eqn:Stein}
\begin{aligned}

X &= {(\mathcal{A}^T-\bar{\mu_i} \mathcal{E}^T)(\mathcal{A}^T+\mu_i \mathcal{E}^T)^{-1}} X_i {(\mathcal{A}^T-\bar{\mu_i} \mathcal{E}^T)^T (\mathcal{A}^T+\mu_i \mathcal{E}^T)^{-T}} \\ &- 2\mu_i (\mathcal{A}^T+\mu_i \mathcal{E}^T)^{-T} \mathcal{C}^T \mathcal{C} (\mathcal{A}^T+\mu_i \mathcal{E}^T)^{-1}

\end{aligned}
\end{equation}

Using the Stein's equation (\ref{eqn:Stein}) in Lemma-\ref{Lres1}, we can find some initial iterates $X_i^{(0)} = 0$ for which the residuals $\mathcal{R}_i^{(k+1)} \ge 0$. Then according to the Lemma-\ref{Lres2} the residual (\ref{eqn:residue2}) can be defined as

\small{
\begin{equation} \label{eqn:residue3}
\begin{aligned}

\mathcal{R}(X_i)& =\left( \prod_{i=1}^{m} {(\mathcal{A}^T-\bar{\mu_i} \mathcal{E}^T)(\mathcal{A}^T+\mu_i \mathcal{E}^T)^{-1}} \right)\mathcal{C}^T
\mathcal{C}\left( \prod_{i=1}^{m} {(\mathcal{A}^T-\bar{\mu_i} \mathcal{E}^T)(\mathcal{A}^T+\mu_i \mathcal{E}^T)^{-1}} \right)^T,\\
&=\mathcal{W}_i \mathcal{W}_i^T,

\end{aligned}
\end{equation}
}

with $\mathcal{W}_i=\left( \prod_{i=1}^{m} {(\mathcal{A}^T-\bar{\mu_i} \mathcal{E}^T)(\mathcal{A}^T+\mu_i \mathcal{E}^T)^{-1}} \right)\mathcal{C}^T$.
\end{proof}

Thus, for $\mathcal{W}_i \mathcal{W}_i^T \leq \tau$ the ADI iterations in the LRCF-ADI algorithm needs to be stopped, where $\tau$ is a given margin of tolerance. Then with the residual factor relation for $V_i$ can be derived as

\begin{equation} \label{eqn:relation1}
\begin{aligned}

V_i =(\mathcal{A}^T+\mu_i \mathcal{E}^T)^{-1} \mathcal{W}_{i-1}.

\end{aligned}
\end{equation}

Then, using (\ref{eqn:relation1}) the residual factor $W_i$ can be derived in a recursive form as

\begin{equation} \label{eqn:relation2}
\begin{aligned}

\mathcal{W}_i=\mathcal{W}_{i-1}-2\Rel(\mu_i) \mathcal{E}^T V_i.

\end{aligned}
\end{equation}

In case of real setting, $\mu_{i+1}:=\bar{\mu_i}$ needs to be considered to find the following form
\begin{equation} \label{eqn:relation3}
\begin{aligned}

\mathcal{W}_{i+1}=\mathcal{W}_{i-1}+2\gamma_i^2 \mathcal{E}^T \left[ \Rel(V_i)+\delta_i \Imag(V_i) \right].

\end{aligned}
\end{equation}

The summary of above techniques is given in the Algorithm-(\ref{pre:LRCF-ADI}).

\begin{algorithm}[t]
\SetAlgoLined
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\caption{Updated LRCF-ADI (Real Version).}
\label{pre:LRCF-ADI}
\Input{ $\mathcal{E}, \mathcal{A}, \mathcal{C}$, $\tau$ (tolerance), $ i_{max}$ (number of iterations) and shift parameters $\{\mu_j\}_{j=1}^{i_{max}}$.}
\Output{Low-rank Cholesky-factor $Z$ such that $X\approx ZZ^T$.}

Consider $\mathcal{W}_0=\mathcal{C}^T, \  Z_0=[ \ ]$ and $i=1$.\\
\While{ $\|\mathcal{W}_{i-1} \mathcal{W}_{i-1}^T\| \geq \tau$ or $i \leq  i_{max}$}{
Solve $V_i=(\mathcal{A}^T+\mu_i \mathcal{E}^T)^{-1} \mathcal{W}_{i-1}$.\\

\eIf{$\Imag(\mu_i)=0$}
{ Update $Z_i=\begin{bmatrix} Z_{i-1} \quad \sqrt{-2\mu_i} V_i \end{bmatrix}$,\\
Compute $\mathcal{W}_i=\mathcal{W}_{i-1}-2\mu_i \mathcal{E}^T V_i$.}
{Assume $\gamma_i=\sqrt{-2\Rel(\mu_i)}$, \quad $\delta_i=\frac{\Rel(\mu_i)}{\Imag(\mu_i)}$,\\
Update $Z_{i+1}=\begin{bmatrix}Z_{i-1} \quad \gamma_i (\Rel(V_i)+\delta_i \Imag(V_i)) \quad \gamma_i \sqrt{\delta_i^2 + 1} \Imag(V_i)\end{bmatrix}$,\\
Compute $\mathcal{W}_{i+1}=\mathcal{W}_{i-1}+2\gamma_i^2 \mathcal{E}^T \left[ \Rel(V_i)+\delta_i \Imag(V_i) \right]$.\\
$i=i+1$
}
$i=i+1$
}
\end{algorithm}

\subsubsection{Adjustment for the unstable systems}

For the unstable index-1 descriptor system, initial feedback matrix $K_0$ needs to be introduced and instead of $(\tilde{\mathcal{A}}^{(i)}, \mathcal{E})$, corresponding shift parameters are needed to be computed from eigen-pair $(\tilde{\mathcal{A}}^{(i)} - \mathcal{B}K_0, \mathcal{E})$. The sparse form of the eigen-pair can be structured as

\begin{equation} \label{eqn:n-eig}
\begin{aligned}

& (\tilde{\mathcal{A}}^{(i)} - \mathcal{B}K_0, \mathcal{E}) =((\mathcal{A}-\mathcal{B}K_0)-\mathcal{B}\mathcal{B}^T (Z^{(i)} (Z^{(i)})^T) \mathcal{E}, \mathcal{E}),\\
&=\left(\begin{bmatrix}  (J_1 - B_1 K_0) - B_1 B_1^T (Z^{(i)} (Z^{(i)})^T) E_1  &  \quad J_2 \\  ( J_3 -  B_2 K_0) - B_2 B_1^T (Z^{(i)} (Z^{(i)})^T) E_1  & \quad J_4 \end{bmatrix}, \begin{bmatrix}  E_1  &  0 \\   0  & 0\end{bmatrix}\right).

\end{aligned}
\end{equation}

To find $V_j^{(i)}$, in each ADI (inner) iteration a shifted linear system needs to be solved as
\begin{equation} \label{eqn:n-linear-st}
\begin{aligned}

((\tilde{\mathcal{A}}^{(i)} -  \mathcal{B}K_0) +\mu_j^{(i)} \mathcal{E})^TV_j^{(i)} &=\mathcal{W}_{j-1}^{(i)},\\
or,((\mathcal{A} - \mathcal{B}K_0) - \mathcal{B}\mathcal{B}^T (Z_{j-1}^{(i)} (Z_{j-1}^{(i)})^T) \mathcal{E} +\mu_j^{(i)} \mathcal{E})^TV_j^{(i)} &=\mathcal{W}_{j-1}^{(i)}.

\end{aligned}
\end{equation}

Thus, $V_j^{(i)}$ can be obtained from the sparse form of the shifted linear system structured as

\small{\begin{equation} \label{eqn:n-broken-st}
\begin{aligned}

\begin{bmatrix}  (J_1  - B_1 K_0 ) - B_1 B_1^T (Z_{j-1}^{(i)} (Z_{j-1}^{(i)})^T) E_1 + \mu_j E_1  &  \quad J_2 \\   (J_3  -  B_2 K_0) - B_2 B_1^T (Z_{j-1}^{(i)} (Z_{j-1}^{(i)})^T) E_1  & \quad J_4 \end{bmatrix}^T \begin{bmatrix}  V_j^{(i)} \\  \Gamma \end{bmatrix}\\ =\begin{bmatrix}  C_1^T  &  E_1^T (Z_{j-1}^{(i)} & (Z_{j-1}^{(i)})^T) B_1 \\  C_2^T   & 0\end{bmatrix}.

\end{aligned}
\end{equation}}

\subsubsection{Estimation of the optimal feedback matrix}

The feedback matrix $K^{(i)}=\mathcal{B}^T X^{(i)} \mathcal{E}=\mathcal{B}^T (Z^{(i)}(Z^{(i)})^T) \mathcal{E}$ needs to be computed in each ADI (inner) iteration and the optimal feedback matrix $K^o=K^{(i_{max})}$ needs to be stored after the final Newton (outer) iteration. The summary of the modified method is given in the Algorithm-(\ref{pre:KN-LRCF-ADI.}).

\begin{algorithm}[]
\SetAlgoLined
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\caption{Modified KN-LRCF-ADI (sparsity-preserving).}
\label{pre:KN-LRCF-ADI.}
\Input{ $E_1, J_1, J_2, J_3, J_4, B_1, B_2, C_1, C_2, K_0$ (initial feedback matrix), and $\tau$ (tolerance).}
\Output{Low-rank Cholesky-factor $Z$ such that $X \approx ZZ^T$ and optimal feedback matrix $K^o$.}

\For {$i \gets 1$ to $i_{max}$}{
Choose $Z_0^{(i)}=[ \  ], $ $K_0^{(i)}=K_0$ and $j=0$.\\
Assume $\mathcal{W}_0^{(i)}= \begin{bmatrix}  C_1^T &  (K^{(i-1)})^T \\ C_2^T & 0 \end{bmatrix}$.\\
Compute adaptive shifts $\left \{  \mu_1^{(i)},......,\mu_J^{(i)} \right \}$ from the eigenpair defined in (\ref{eqn:n-eig}).\\
\While{ $\left(\|\mathcal{W}_j^{(i)}\|^2> \tau \|\mathcal{W}_0^{(i)}\|^2 \right)$}{
$j=j+1$\\
Solve the linear system (\ref{eqn:n-broken-st}) for $V_j^{(i)}$.\\

\eIf{$\text{Im}(\mu_j^{(i)})=0$}
{Update $Z_j^{(i)}=\begin{bmatrix} Z_{j-1}^{(i)} \quad & \sqrt{-2 \mu_i} V_j^{(i)} \end{bmatrix}$,\\
Compute $\mathcal{W}_j^{(i)}=\mathcal{W}_{j-1}^{(i)}-2 \mu_j^{(i)} E_1^T V_j^{(i)}$,\\
Compute $K_j^{(i)}=K_{j-1}^{(i)}-2 \mu_j^{(i)} (B_1 - J_2{J_4}^{-1}B_2)^T V_j^{(i)} (V_j^{(i)})^T E_1$.}
{Assume $\gamma_j^{(i)}=\sqrt{-2\text{Re}(\mu_j^{(i)})}$, \quad $\beta_j^{(i)}=\frac{\text{Re}(\mu_j^{(i)})}{\text{Im}(\mu_j^{(i)})}$,
\quad $\delta_j^{(i)} = \text{Re}(V_j^{(i)}) + \beta_j^{(i)}\text{Im}(V_j^{(i)}),$\\
Compute $ Z_d=\begin{bmatrix} \gamma_j^{(i)} \delta_j^{(i)} \quad & \gamma_j^{(i)} \sqrt{(\beta_j^{(i)})^2+1} \ \text{Im}(\mu_j^{(i)}) \end{bmatrix}$,\\
Update $Z_{j+1}^{(i)}=\begin{bmatrix} Z_{j-1}^{(i)}  \quad & Z_d \end{bmatrix}$,\\
Compute $\mathcal{W}_{i+1}^{(i)} = \mathcal{W}_{i-1}^{(i)} -4\text{Re}(\mu_j^{(i)}) E_1^T \delta_j^{(i)}$,\\
Compute $K_{j+1}^{(i)}=K_{j-1}^{(i)}+(B_1 - J_2{J_4}^{-1}B_2)^T Z_d {Z_d}^T E_1$,\\
$j=j+1$.
}
}
Update $Z^{(i)}=Z_j^{(i)}$ and $K^{(i)}=K_j^{(i)}$.
}
\end{algorithm}

\subsection{Computation of the optimally stabilized system and the optimal control}

The optimal feedback matrix $K^o=\mathcal{B}^T X\mathcal{E}$ can be achieved by the feasible solution $X$ of the Riccati equation (\ref{eqn:D-GCARE}). Then, applying $\mathcal{A}_s = \mathcal{A} - \mathcal{B} K^o$, optimally stabilized LTI continuous-time system can be written as (\ref{eqn:st-state-space}). To preserve the structure of the system, it needs to back to the oginial structure (\ref{eqn:DS_matrix_vector}), and for this the submatrices $J_1, J_3$ are replaced by $J_1 - B_1 K^o, J_3 - B_2 K^o$, respectively.

Finally, the desired optimal control $u^o = -K^o x(t)$ of the targeted power system models can be computed.

\section{Numerical Results}

The stability of the target models is investigated and the unstable models are stabilized through the Riccati-based feedback stabilization process. The proposed methods are employed to find the solution of the Riccati equation arising from the BIPS models and corresponding feedback matrices are generated for system stabilization. Also, initial Bernoulli feedback stabilization is implemented for the convenient rate of convergence.
All the results have been achieved using the MATLAB 8.5.0 (R2015a) on a Windows machine having Intel-Xeon Silver 4114 CPU \@$2.20$ GHz clock speed, $2$ cores each and $64$ GB of total RAM.

\subsection{Brazilian inter-connected power system models}\label{bips}

Power system models are an essential part of engineering fields that consists of simulations based on power generations and grid networks. The computation required to analyze electrical power systems employing mathematical models utilizing real-time data. There are several applications of the power system model, i.e., electric power generation, utility transmission and distribution, railway power systems, and industrial power generation \cite{martins2007computation}.

The Brazilian Inter-connected Power Systems (BIPS) is one of the most convenient examples of the power system models with various test systems \cite{freitas1999computationally}. The following Table~\ref{tab:BIPS_systems} provides the details about the models. The detailed structure of those will be found at \textit{https://sites.google.com/site/rommes/software}, where all of them are index-1 descriptor system. The models $mod-606$, $mod-1998$, $mod-2476$ and $mod-3078$ have the unstable eigenvalues, whereas the models $mod-1142$, $mod-1450$ and $mod-1963$ have stable eigenvalues \cite{leandro2015identification}. Here the names of the models are considered according to their number of states.

\begin{table}[H]
\centering
\caption{Structure of the Models derived from BIPS test systems}
\label{tab:BIPS_systems}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
\begin{tabular}[c]{@{}c@{}}Test \\ systems\end{tabular} & Dimensions & States & \begin{tabular}[c]{@{}c@{}}Algebraic\\ variables\end{tabular} & Inputs & Outputs \\ \hline
\multirow{4}{*}{BIPS98} & 7135 & 606 & 6529 & 4 & 4 \\ \cline{2-6}
		& 9735 & 1142 & 8593 & 4 & 4 \\ \cline{2-6}
		& 11265 & 1450 & 11582 & 4 & 4 \\ \cline{2-6}
		& 13545 & 1963 & 13068 & 4 & 4 \\ \hline
\multirow{3}{*}{BIPS07} & 15066 & 1998 & 13068 & 4 & 4 \\ \cline{2-6}
		& 16861 & 2476 & 14385 & 4 & 4 \\ \cline{2-6}
		& 21128 & 3078 & 18050 & 4 & 4 \\ \hline
\end{tabular}
\end{table}

\subsection{Comparison of the results found by RKSM and KN-LRCF-ADI}

The CAREs arising from the models $mod-606$, $mod-1998$ and $mod-2476$ are efficiently solved and stabilized the corresponding models by both RKSM and KN-LRCF-ADI techniques. As model $mod-3078$ is semi-stable, the computation of CARE derived from this model is not possible by LRCF-ADI techniques but by the RKSM approach model $mod-3078$ successfully stabilized and the numerical result for model $mod-3078$ is investigated for RKSM only.

\begin{table}[H]
\centering
\caption{Results of RKSM applied BIPS models}
\label{tab:rksm}
\begin{tabular}{|c|c|c|c|c|}\hline
Model & \begin{tabular}[c]{@{}c@{}}No of \\ iterations\end{tabular} & Tolerance & \begin{tabular}[c]{@{}c@{}}Numerical\\ rank\end{tabular} & \begin{tabular}[c]{@{}c@{}}CPU time\\ (second)\end{tabular} \\ \hline
		$mod-606$ & 100 & $10^{-10}$ & 195 & $1.81 \times 10^2$ \\ \hline
		$mod-1998$ & 200 & $10^{-10}$ & 266 & $1.41 \times 10^3$ \\ \hline
		$mod-2476$ & 248 & $10^{-10}$ & 265 & $3.06 \times 10^3$ \\ \hline	$mod-3078$ & 257 & $10^{-5}$ & 295 & $3.01 \times 10^3$ \\ \hline
\end{tabular}
\end{table}

Table~\ref{tab:rksm} depicts the numerical results of the stabilization process via RKSM for the unstable BIPS models and various properties of the stabilized systems are illustrated, whereas Table~\ref{tab:kn-adi} displays the several modes of ADI techniques in KN-LRCF-ADI method for stabilizing the unstable BIPS models including characteristics of the stabilized models. In both of the tables, Table~\ref{tab:rksm} and Table~\ref{tab:kn-adi} we have analyzed the same features of the stabilized BIPS models, we can easily compare the efficiency and robustness of the proposed methods.

From the above tables it can be said that the proposed RKSM approach has quick convergence ability and occupies very small solution space to provide the efficient solution of the CAREs. In contrast LRCF-ADI based Kleinman-Newton has several approaches for finding the solution of CAREs, whereas most of the approaches required higher computation time. Also, there are deviations in the numerical ranks of the factored solution of CAREs in the Kleinman-Newton approaches and in all of the cases RKSM provides significantly better result.

\begin{table}[H]
\centering
\caption{Results of KN-LRCF-ADI applied BIPS models}
\label{tab:kn-adi}
\begin{tabular}{|c|c|c|c|c|}\hline
Model & \begin{tabular}[c]{@{}c@{}}Total\\ iterations\end{tabular} & Tolerance & \begin{tabular}[c]{@{}c@{}}Numerical\\ rank\end{tabular} & \begin{tabular}[c]{@{}c@{}}CPU time\\ (second)\end{tabular} \\ \hline
\multirow{4}{*}{$mod-606$} & 311 & $10^{-5}$ & 481 & $2.74 \times 10^2$ \\ \cline{2-5}
		& 544 & $10^{-10}$ & 953 & $7.21 \times 10^2$ \\ \cline{2-5}
		& 508 & $10^{-5}$ & 473 & $4.10 \times 10^2$ \\ \cline{2-5}
		& 853 & $10^{-10}$ & 969 & $7.92 \times 10^2$ \\ \hline
\multirow{4}{*}{$mod-1998$} & 277 & $10^{-5}$ & 663 & $2.30 \times 10^3$ \\ \cline{2-5}
		& 485 & $10^{-10}$ & 1201 & $5.83 \times 10^3$ \\ \cline{2-5}
		& 514 & $10^{-5}$ & 497 & $3.27 \times 10^3$ \\ \cline{2-5}
		& 1003 & $10^{-10}$ & 1417 & $1.17 \times 10^4$ \\ \hline
\multirow{4}{*}{$mod-2476$} & 254 & $10^{-5}$ & 473 & $3.13 \times 10^3$ \\ \cline{2-5}
		& 363 & $10^{-10}$ & 937 & $5.75 \times 10^3$ \\ \cline{2-5}
		& 366 & $10^{-5}$ & 441 & $3.02 \times 10^3$ \\ \cline{2-5}
		& 849 & $10^{-10}$ & 849 & $9.50 \times 10^4$ \\ \hline
\end{tabular}
\end{table}

\subsection{Stabilization of eigenvalues}

Fig.\ref{fig:eig_606}, Fig.\ref{fig:eig_1998}, and Fig.\ref{fig:eig_2476} illustrates the eigenvalue stabilization of the models $mod-606$, $mod-1998$, $mod-2476$ by both the RKSM and LRCF-ADI based Kleinamn-Newton techniques.

\begin{figure}[H]
	\centering
	\begin{subfigure}{.95\linewidth}
		\includegraphics[width=\linewidth]{figures1/uddin1.PNG}
		\vspace{0.01in}
		\caption{Stabilized by RKSM}
		\vspace{0.1in}
		\label{fig:eig606_1}
	\end{subfigure}
	\vspace{0.1in}
	\begin{subfigure}{.95\linewidth}
		\includegraphics[width=\linewidth]{figures1/uddin2.PNG}
		\vspace{0.01in}
		\caption{Stabilized by KN-LRCF-ADI}
		\vspace{0.1in}
		\label{fig:eig606_2}
	\end{subfigure}
	\vspace{0.1in}
	\caption{Comparisons of the eigenvalues for the model $mod-606$}
	\label{fig:eig_606}
\end{figure}

\begin{figure}[H]
	\centering
	\begin{subfigure}{.95\linewidth}
		\includegraphics[width=\linewidth]{figures1/uddin3.PNG}
		\vspace{0.01in}
		\caption{Stabilized by RKSM}
		\vspace{0.1in}
		\label{fig:eig1998_1}
	\end{subfigure}
	\vspace{0.1in}
	\begin{subfigure}{.95\linewidth}
		\includegraphics[width=\linewidth]{figures1/uddin4.PNG}
		\vspace{0.01in}
		\caption{Stabilized by KN-LRCF-ADI}
		\vspace{0.1in}
		\label{fig:eig1998_2}
	\end{subfigure}
	\vspace{0.1in}
	\caption{Comparisons of the eigenvalues for the model $mod-1998$}
	\label{fig:eig_1998}
\end{figure}

From the sub-figures in the above mentioned figures, it has been observed that the eigenvalues of the models $mod-606$ and $mod-2476$ are stabilized very well but its marginal for the model $mod-1998$. Thus, it can be concluded that both the RKSM and LRCF-ADI based Kleinamn-Newton techniques have adequate efficiency to stabilize the unstable descriptor systems by closed-loop structure via Riccati based feedback stabilization.

Fig.\ref{fig:eig_3078} displays the applicability of the RKSM for the semi-stable descriptor system, whereas LRCF-ADI based methods are ineffective in this case. Since the eigenvalues of the semi-stable model $mod-3078$ are very close to the imaginary axis, it is not possible to realize in the normal view. So, for the simulation tool capacity and visual convenience, we have considered the a very magnified view of the eigenspaces.

\begin{figure}[H]
	\centering
	\begin{subfigure}{.95\linewidth}
		\includegraphics[width=\linewidth]{figures1/uddin5.PNG}
		\vspace{0.01in}
		\caption{Stabilized by RKSM}
		\vspace{0.1in}
		\label{fig:eig2476_1}
	\end{subfigure}
	\vspace{0.1in}
	\begin{subfigure}{.95\linewidth}
		\includegraphics[width=\linewidth]{figures1/uddin6.PNG}
		\vspace{0.01in}
		\caption{Stabilized by KN-LRCF-ADI}
		\vspace{0.1in}
		\label{fig:eig2476_2}
	\end{subfigure}
	\vspace{0.1in}
	\caption{Comparisons of the eigenvalues for the model $mod-2476$}
	\label{fig:eig_2476}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=\linewidth]{figures1/uddin7.PNG}
	\vspace{0.01in}
	\caption{Comparisons of the eigenvalues (stabilized by RKSM) for the model $mod-3078$}
	\vspace{0.1in}
	\label{fig:eig_3078}
\end{figure}

\subsection{Stabilization of step-responses}

The investigation of the figures from Fig.\ref{fig:step-res_606_1}-Fig.\ref{fig:step-res_3078_2} consist the step-responses for some dominant input/output relations to compare the RKSM and the LRCF-ADI based Kleinman-Newton approaches via the system stabilization. Since the power system models are of $4\times4$ input-output relations, there are $16$ step-responses for each models. For the effective comparison we have investigated only some graphically significant step-responses.

\begin{figure}[H]
	\centering
	\includegraphics[width=\linewidth]{figures1/uddin8.PNG}
	\vspace{0.01in}
	\caption{Comparisons of step-responses for the model $mod-606$ for the first-input/third-output}
	\label{fig:step-res_606_1}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=\linewidth]{figures1/uddin9.PNG}
	\vspace{0.01in}
	\caption{Comparisons of step-responses for the model $mod-606$ for the second-input/first-output}
	\label{fig:step-res_606_2}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=\linewidth]{figures1/uddin10.PNG}
	\vspace{0.01in}
	\caption{Comparisons of step-responses for the model $mod-1998$ for the third-input/second-output}
	\label{fig:step-res_1998_1}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=\linewidth]{figures1/uddin11.PNG}
	\vspace{0.01in}
	\caption{Comparisons of step-responses for the model $mod-1998$ for the fourth-input/third-output}
	\label{fig:step-res_1998_2}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=\linewidth]{figures1/uddin12.PNG}
	\vspace{0.01in}
	\caption{Comparisons of step-responses for the model $mod-2476$ for the first-input/third-output}
	\label{fig:step-res_2476_1}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=\linewidth]{figures1/uddin13.PNG}
	\vspace{0.01in}
	\caption{Comparisons of step-responses for the model $mod-2476$ for the third-input/second-output}
	\label{fig:step-res_2476_2}
\end{figure}

From Fig.\ref{fig:step-res_606_1}-Fig.\ref{fig:step-res_2476_2} it is evident that both RKSM and LRCF-ADI integrated Kleinman-Newton techniques are applicable for the Riccati based feedback stabilization of unstable power system models. The graphical comparisons indicates by RKSM is suitably robust. On the other hand, though sometimes the Kleinman-Newton approach provides very good accuracy but it has some scattered behaviors.

Moreover, Fig.\ref{fig:step-res_3078_1} and Fig.\ref{fig:step-res_3078_2} shows the applicability of the RKSM technique for the Riccati based feedback stabilization for the semi-stable index-1 descriptor systems.

\begin{figure}[H]
	\centering
	\includegraphics[width=\linewidth]{figures1/uddin14.PNG}
	\vspace{0.1in}
	\caption{Comparisons of step-responses for the model $mod-3078$ for the second-input/first-output}
	\label{fig:step-res_3078_1}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=\linewidth]{figures1/uddin15.PNG}
	\vspace{0.1in}
	\caption{Comparisons of step-responses for the model $mod-3078$ for the fourth-input/third-output}
	\label{fig:step-res_3078_2}
\end{figure}

\section{Conclusion}

From the tabular and graphical comparisons of the results of numerical computations, we have observed that by both RKSM and KN-LRCF-ADI techniques, the CAREs arising from the unstable index-1 descriptor systems are efficiently solved and the corresponding models are properly stabilized. The semi-stable index-1 descriptor system is successfully stabilized through Riccati-based feedback stabilization by RKSM, whereas KN-LRCF-ADI is still not suitable for it. There are deviations of the numerical ranks of the factored solutions of CAREs in the Kleinman-Newton approaches, while RKSM provides significantly better results for all the cases. RKSM approach has quick convergence ability and occupies very small solution spaces to provide efficient solutions to the CAREs. In contrast, LRCF-ADI based Kleinman-Newton has several approaches for finding the solutions of CAREs, where almost all of the approaches required higher computation time. Riccati-based feedback stabilization for the index-1 descriptor systems by the RKSM approach is very effective and robust. Contrariwise, LRCF-ADI based Kleinman-Newton method is slightly scattered in case of the stabilization of step-responses. Thus, it can be concluded that the RKSM is suitably applicable to the unstable index-1 descriptor systems for Riccati-based feedback stabilization and this method is more preferable to the Kleinman-Newton method in the sense of computation time and memory allocation.

In this work, the MATLAB library command \texttt{care} is used in RKSM to find the solution of the CAREs governed from the reduced-order models and used the direct backward inversion technique for solving shifted sparse linear systems. In the future, we will try to find the self-sufficient RKSM algorithm for solving CAREs and apply the "Restarted Hessenberg Method" for solving shifted sparse linear systems.  Since the LRCF-ADI techniques are inefficient for the semi-stable systems, we will work for it as well.

\section{Acknowledgment}

An earlier version of it has been presented as a pre-print in "Riccati-based feedback stabilization for unstable power system models". The earlier version of this work can be found at the link https://arxiv.org/abs/2006.14210.

This work is under the project \textit{"Computation of Optimal Control for Differential-Algebraic Equations (DAE) with Engineering Applications"}. This project is funded by United International University, Dhaka, Bangladesh. It starts from October 01, 2019, and the reference is IAR/01/19/SE/18.

\bibliographystyle{IEEEtran}
\bibliography{Ref_MUn}

\end{document}

config: arxiv


text:
\section*{Acknowledgements}
\label{sec:ack}

This work was partially supported by
ANPCyT PICT-2015-2218,
UBACyT Grants
	20020190100126BA,
	20020170100495BA and
	20020160100095BA
(Argentina),
and
Programa Regional MATHAMSUD MATH190013.
Carolina L. Gonzalez is partially supported by a CONICET doctoral fellowship.


TITLE: Given a specific function $f$, how to compute the left-inverse of $f$ in the sense of $\approx$?
QUESTION [2 upvotes]: For a non-negative function $\varphi$ defined on $[0,\infty)$, the left-inverse $\varphi^{-1}$ of $\varphi$ is defined by setting, $\forall t\geq 0$,
$$\varphi^{-1}(t):=\inf\{u\geq0:\varphi(u)\geq t\}.$$
For two functions $\varphi$ and $\psi$, we say that $\varphi\approx\psi$, if there exist constants $c_1,c_2>0$, such that $\forall t\geq0$,
$$c_1\varphi(t)\leq\psi(t)\leq c_2\varphi(t).$$
Let $p\geq1$, define $f(t):=t^p\ln(e+t)$, $t\geq0$ (this is an example appearing in the paper about equations). How to compute the $f^{-1}$ in the sense of $\approx$? I guess its
answer is $\forall t\geq0$, $$f^{-1}(t)\approx\left[\frac{t}{\ln(e+t)}\right]^{1/p}\quad \text{or}\quad f^{-1}(t)\approx \min\left\{t^{1/p},\left[\frac{t}{\ln(e+t)}\right]^{1/p}\right\},$$ but I do not know how to compute $f^{-1}$ in the sense of $\approx$, I need some help.

REPLY [0 votes]: As a supplement and extension of Iosif's answer, let me mention that, more generally, if $$f(t) \sim t^{pq} (\ell(t^q))^p \qquad \text{as } t \to \infty$$ for a slowly varying function $\ell$, then $$f^{-1}(s) \sim s^{1/(pq)} (\ell^\#(s^{1/p}))^{1/q} \qquad \text{as } s \to \infty,$$ where $\ell^\#$ is another slowly varying function: the de Bruijn conjugate of $\ell$. See Proposition 1.5.15 in Regular Variation by Bingham, Goldie and Teugels. If $\ell(t) = \log t$, then $\ell^\#(s) \sim 1 / \log t$, and so we recover an asymptotic variant of the desired result (with $q = 1$).
A variant with asymptotic equality ($\sim$) replaced by comparability ($\approx$) is essentially given as Problem 1.11.19.


TITLE: Existence of imaginary quadratic fields of class numbers coprime to $p$ with prescribed splitting behaviour of $p$
QUESTION [7 upvotes]: Let $x\in\{\text{totally ramified, inert, totally split}\}.$
If $p\geq 5$ is a prime, are there infinitely many imaginary quadratic fields $K=\mathbb{Q}(\sqrt{-d})$ of class number coprime to $p$ so that $p$ has ramification behaviour $x$ in $K/\mathbb{Q}$?

REPLY [10 votes]: A more general version of this statement was shown by Kimura in Acta Arith. (2003).  His corollary gives $\gg \sqrt{X}/\log X$ such fields ${\Bbb Q}(\sqrt{-d})$ with $d\le X$, and also allows you to add further
splitting conditions.  There is an extensive literature on divisibility and indivisibility of class numbers and Kimura's paper has many other relevant references.
Update  There is also an erratum to Kimura's paper -- it seems that he needs the existence of one such field to get a lower bound for the number of such fields; for large $p$ this is guaranteed by work of Horie.  In the meantime a paper of Wiles has addressed the general version of this question, and there is a quantification of Wiles's result due to Beckwith.


TITLE: Prove that the following are equivalent.
QUESTION [0 upvotes]: Let $(V,d)$ be a metric space and $A \subset V$.  Define $f \colon V \to \mathbb{R}$ by $f(x) = 1$ if $x\in A$ and $f(x) = 0$ if $x\notin A$. Prove the following are equivalent:

$A$ is open and closed.
$f$ is continuous.

I have the feeling that i have to prove that $f(A)$ is also open and closed but I am not quite sure.

REPLY [2 votes]: I think you mean that  $f(x) = 0$ if $x\notin A$.
In that case you probably have to use the definition of $f continuous$ in a metric set. $f$ is continuous iff for all open (resp closed) set $B$, $f^{-1}(B)$ is open (resp closed).
$(1) \implies (2)$:
If $A$ is open and closed, then let $B$ an open set of $R$.
If $0, 1 \in B$ then $f^{-1}(B) = V$ which is open.
If $0, 1 \notin B$ then $f^{-1}(B) = \varnothing$ which is open.
If $0\in B$ and $1 \notin B$ then $f^{-1}(B) = V\backslash A$
which is open because A is closed.
If $1\in B$ and $0 \notin B$ then f^{-1)(B) = A$ which is open .
So $f$ is continuous.
$(2) \implies (1)$:
If $f$ is continuous, then $A = f^{-1}(\{1\})$ is closed because ${1}$ is closed in $\mathbb{R}$. Moreover,  $A = f^{-1}(]-1/2;1/2[)$ is open because $]-1/2;1/2[$ is open in $\mathbb{R}$

REPLY [0 votes]: hint for $2)\implies 1) $
$f $ continuous and $\{1\} $ is closed $\implies f^{-1}(\{1\}) =A$ is closed.
$f $ continuous and $\{0\} $ closed $\implies f^{-1}(\{0\})=V\backslash A $ closed $\implies A $open.


TITLE: having trouble with limit of integral
QUESTION [7 upvotes]: How do I solve the following?
$$\lim_{x\to 0} \int_0^1 \cos\left(\frac{1}{xt}\right)\, dt$$

REPLY [4 votes]: Here's another take.  Since $\cos \left(\frac{1}{xt}\right)$ is an even function of $x$, we can take the limit from the right or from the left, and the result is the same.  So assume $x > 0$.  Then apply the $u = 1/(xt)$ substitution as in joriki's answer, followed by integration by parts in the opposite direction from joriki's answer.  Rewriting the result in terms of the sine integral produces the output from Wolfram Alpha mentioned by Eivind and Sivaram:
$$ \lim_{x \to 0} \left(\frac{1}{x} \int_0^{1/x} \frac{\sin u}{u} du  - \frac{\pi}{2x} + \cos \left(\frac{1}{x}\right)\right).$$
Now, apply the following known series expansion for the sine integral $\int_0^x \frac{\sin u}{u} du$  (valid for large values of $x$):
$$\int_0^x \frac{\sin u}{u} du = \frac{\pi}{2} -\frac{\cos x}{x} \left(1 - \frac{2!}{x^2} + \frac{4!}{x^4} \pm \cdots \right) - \frac{\sin x}{x} \left(\frac{1}{x} - \frac{3!}{x^3} \pm \cdots \right).$$
After simplifying, we are left with
$$\lim_{x \to 0} \left(-\cos \left(\frac{1}{x}\right) \left(- 2! x^2 + 4! x^4 \pm \cdots \right) - \sin\left(\frac{1}{x}\right) \left(x - 3! x^3 \pm \cdots \right) \right).$$
Since $\cos \left(\frac{1}{x}\right)$ and $\sin \left(\frac{1}{x}\right)$ are both bounded by $-1$ and $1$, the limit is $0$.
In addition, we can see that the dominant term in the limit is $- x\sin\left(\frac{1}{x}\right)$, which is clearly bounded by $|x|$, as conjectured by Sivaram in the comments.


TITLE: Functional differential equation $D'(t) = 2 D(2t) - D(t)$
QUESTION [7 upvotes]: Can the functional differential equation $$D'(t) = 2 D(2t) - D(t)$$
with initial condition $D(0)=0$
be solved in closed form?
A Laplace transform yields the functional equation
$$
(s+1)d(s) = d(s/2)
$$
but this equation seems badly behaved ($d(0)$ indeterminate).

REPLY [6 votes]: There  are series solutions
$$ \eqalign{d(t) &= c \left( e^{-t} - 2 e^{-2t} + \frac{4}{3} e^{-4t} - \frac{8}{21} e^{-8t} + \ldots\right)\cr
                 &= c \sum_{k=0}^\infty \frac{(-2)^k}{ \prod_{j=1}^k (2^j-1)} e^{-2^kt}\cr}$$
which all have $d(0) = 0$. These should converge for $t \ge 0$.

REPLY [3 votes]: Let us relax the initial condition and look at a solution of the form
$$ D(t) = \sum_{k=-\infty}^{\infty} a_k \exp(-2^k t). $$
Plugging this to the equation and assuming that term-wise differentiation applies, we find that $(a_k)$ satisfies the recurrence relation $ (1-2^k) a_k = 2a_{k-1}$, whose solution can be written in terms of the $q$-Pochhammer symbol:
$$ a_k = \begin{cases} \dfrac{2^k}{(2,2)_k} a_0, & k \geq 0 \\ 0, & k < 0 \end{cases} $$
Consequently we can write
$$ D(t) = a \sum_{k=0}^{\infty} \frac{2^k}{(2;2)_k}\exp(-2^k t) $$
for some constant $a$.
p.s. As pointed out by Robert Israel, these solutions indeed satisfy the initial condition $D(0) = 0$ (which I did not expect before looking at his answer). Here is a simple confirmation:
\begin{align*}
D(0)
&= a \sum_{k=0}^{\infty} \frac{2^k}{(2;2)_k} \\
&= a \left(1 + \sum_{k=1}^{\infty} \frac{2^k - 1}{(2;2)_k} + \sum_{k=1}^{\infty} \frac{1}{(2;2)_k} \right) \\
&= a \left(1 - \sum_{k=1}^{\infty} \frac{1}{(2;2)_{k-1}} + \sum_{k=1}^{\infty} \frac{1}{(2;2)_k} \right) \\
&= 0.
\end{align*}


TITLE: Find the singularities of $f(z) =\frac{1}{(2\sin z - 1)^2}$.
QUESTION [3 upvotes]: Find the singularities of $f(z) =\frac{1}{(2\sin z - 1)^2}$.
I am just learning about singularities and I was wondering if someone could give me feedback on my work. So I think, for this function, that there are singularities at $z=\frac{\pi}{6}+2k\pi,\frac{5\pi}{6}+2k\pi$, would this be correct? Additionally, I am classifying them as essential and nonremovable. Is this also correct?
To determine that they were nonremovable, I tool the limit of the function approaching the singularities and found that they tended towards infinity.

REPLY [2 votes]: This $f$ has a singularity anywhere \begin{align*}
2 \sin z - 1 &= 0  \text{,}  \\
\sin z &= 1/2  \text{.}
\end{align*}
This has infinitely many solutions in the reals,
\begin{align*}
z \in &\{ \arcsin(1/2) + 2 \pi k  : k \in \Bbb{Z} \} \cup  \\
\quad &\{\pi - \arcsin(1/2) + 2 \pi k : k \in  \Bbb{Z}\}  \\
 = & \{ \pi/6 + 2 \pi k  : k \in \Bbb{Z} \} \cup  \\
\quad &\{5\pi/6 + 2 \pi k : k \in  \Bbb{Z}\}  \text{.}
\end{align*}
It turns out these are the only two families of solutions, but you have done nothing to show that there are not more solutions.  To do so, start with, for $x,y\in\Bbb{R}$,
$$  \sin(x+\mathrm{i} \, y) = \sin(x) \cosh(y) + \mathrm{i} \cos(x) \sinh(y)  \text{.}  $$
Then talk about the zeroes of the real-valued cosine and hyperbolic sine (since $1/2$ has imaginary part $0$), which will restrict $y$ so that $\cosh = 1$, forcing the real solutions to be the only solutions.
If you have shown that the limits are infinity approaching each of these infinitely many singularities, then you have shown that they are not removable.  I'm not convinced you have understood the definition/nature of essential singularities.