# Constrained velocity-free control of spacecraft attitude via explicit
reference governor

Qingqing Dang  Wenbo Li111<EMAIL_ADDRESS>Haichao Gui

###### Abstract

This paper introduces an explicit reference governor-based control scheme
tailored for addressing the velocity-free spacecraft attitude maneuver
problem. This problem is subject to specific constraints, namely the pointing
constraint, angular velocity constraint, and input saturation. The proposed
control scheme operates in two layers, ensuring the asymptotic stability of
the spacecraft’s attitude while adhering to the aforementioned constraints.
The inner layer employs output feedback control utilizing an angular velocity
observer based on immersion and invariance technology. This observer
facilitates attitude stabilization without the measurement of angular
velocity. Through an analysis of the geometry associated with the pointing
constraint, determination of the upper bound of angular velocity, and
optimization of the control input solution, the reference layer establishes a
safety boundary described by the invariant set. Additionally, we introduce the
dynamic factor related to the angular velocity estimation error into the
invariant set to prevent states from exceeding the constraint set due to
unmeasurable angular velocity information. The shortest guidance path is then
designed in the reference layer. Finally, we substantiate the efficacy of the
proposed constrained attitude control algorithm through numerical simulations.

###### keywords:

Attitude maneuver, Constraints , Explicit reference governor, Angular velocity
observer

††journal:  ArXiv

[label1] organization=School of Civil Aviation, Northwestern Polytechnical
University,city=Xi’an, postcode=710072, state=Shanxi, country=China

[label2]organization=Beijing Institute of Control Engineering, China Academy
of Space Technology,city=Beijing, postcode=100190, country=China

[label3]organization=School of Astronautics, Beihang University, city=Beijing,
postcode=100191, country=China

##  1 Introduction

Spacecraft attitude maneuver plays a significant role in complex space
autonomous missions[1]. Restricted by the actuators and sensitive payloads,
attitude maneuver algorithms must ensure system stability while simultaneously
adhering to multiple constraints [2, 3, 4]. For instance, the spacecraft is
usually required to maneuver from one state to another within the defined time
, while keeping its star sensor avoid from the bright objects (e.g. earth) and
preventing the command torque from exceeding the capacity of the actuator [5].
These missions exemplify spacecraft maneuvering amidst state and control
constraints [6, 7]. Moreover, in scenarios where gyroscopes fail, the
unavailability of angular velocity information becomes a significant challenge
[8]. Consequently, constrained velocity-free attitude control is an issue of
great theoretical and practical importance.

For attitude control systems with actuator saturation, if the input
constraints are not considered in the controller design explicitly, although
the performance is affected by the input limitation, its stability sometimes
can still be proved theoretically [9]. Traditional controllers designed
directly through Lyapunov function lack the ability to restrict state
trajectories. Hence, the attitude commands are used in the attitude maneuver
path design in the presence of multiple attitude constraints [5, 10]. This
strategy can effectively solve part of the engineering problems, but it
suffers from limited flexibility and struggles to meet tasks demanding high
real-time dynamic requirements. Addressing challenges caused by actuator
saturation, control bandwidth limits, slew rate constraints, and/or eigenaxis
slew constraints, Wie et al. introduced saturation and integration functions
within a nonlinear feedback control logic for rapid re-targeting control of
agile spacecraft [11]. This method can handle single-axis maneuvers with
particular constraints well. However, it encounters difficulties in handling
three-axis maneuvers with intricate constraints.

The amalgamation of potential functions and Lyapunov functions presents a
promising technology for addressing complex constraints in the constrained
attitude control problem [12, 13, 14]. Lee et al. [15] constructed a strictly
convex logarithmic barrier potential for attitude-constrained zones by
utilizing a convex parameterization technology. Inspired by[15] and using the
anti-unwinding attitude error function, a new algorithm for the attitude
reorientation guidance under forbidden pointing constraints is proposed in
[16]. Furthermore, Shen et al. [17] addressed the rest-to-rest three-axis
attitude reorientation under multiple attitude-constraint zones and angular
velocity limits via a quadratic potential function and a logarithmic potential
function. Nevertheless, it is difficult to simultaneously handle different
types of complex constraints by the potential functions based constrained
control algorithm. Since the potential function is constructed in the Lyapunov
function and the convergence of Lyapunov function is the result of the
convergence game between potential function and states, the robustness of the
system may become worse.

Trajectory optimization methods, notably model predictive control (MPC), offer
a means to tackle constrained control issues. In [18] and[19], MPC on SO(3)
has been developed for constrained attitude maneuver of fully actuated
spacecraft. However, the necessity to optimize the function at each sampling
horizon in MPC restricts its application in systems requiring rapid response,
such as spacecraft maneuvering. Recently, a novel add-on control scheme called
explicit reference governor (ERG) was introduced by Nicotra et al. [20, 21,
22, 23]. The key idea is to augment a pre-stabilized system with a control
unit and manipulates the auxiliary reference to ensure constraint
satisfaction, which means the stability and the constraint issues can be
handled separately. This control technology has found application in solving
problems related to Unmanned Aerial Vehicles and spacecraft attitude control
with state constraints [3, 7, 24, 25].

In addition to constraints, another challenge in attitude control is the
velocity-free control problem. This issue has drawn significant attention from
researchers and has been extensively studied [26, 27, 28]. For instance,
immersion and invariance (I&I) technology have been utilized to develop a
globally exponentially convergent observer for the angular velocity in [9, 27,
29, 30]. In our earlier research, a six-degree-of-freedom observer was
constructed using the I&I approach [31]. However, the velocity-free attitude
maneuver problem in the presence of constraints was studied in just a few
works. For instance, a velocity-free attitude reorientation control law with
pointing constraints is established in [32].

Inspired by the ERG and the I&I technologies, a constrained velocity-free
control algorithm for spacecraft reorientation is presented in this paper. The
algorithm takes into account attitude pointing, angular velocity, and control
input constraints. The formulation of attitude dynamics and various
constraints is expressed using modified Rodrigues parameters (MRPs). The MRPs
constitute a singular, nonunique and minimal parametrization set of the three-
dimensional special orthogonal group SO(3). Fortunately, the singularity can
be avoided by using the nonuniqueness properties through switching the
parameters between MRPs and its shadow at the unit sphere [33, 34].
Subsequently, the ERG-based control scheme is derived, wherein the output
controller, relying on the angular velocity observer, is primarily designed in
the inner loop. To the best of the authors’ knowledge, the result presented in
this paper is the first attempt to address the observer-based attitude
maneuver issue with pointing constraints, angular velocity constraints, and
input constraints. Finally, the performance and robustness of the proposed
algorithm is verified by the numerical simulations and Monte Carlo
simulations.

##  2 Preliminaries

###  2.1 Spacecraft attitude kinematics and dynamics

The MRPs vector is defined in terms of an Euler rotation angle
$\phi\in\mathbb{R}$ about the principal axis
$\\{\textbf{{n}}|\textbf{{n}}^{T}\textbf{{n}}=1,\textbf{{n}}\in\mathbb{R}^{3}\\}$.
Let $\mathscr{F}_{B}$ be the body-fixed frame, and $\mathscr{F}_{I}$ be the
inertial frame. Then, the attitude with respect to the inertial frame can be
described by MPRs and given by
$\boldsymbol{\sigma}_{BI}=\textbf{{n}}_{BI}{\rm{tan}}(\phi_{BI}/4)$. The
attitude kinematics and dynamics of the rigid-body spacecraft are given by
[35]

$\dot{\boldsymbol{\sigma}}_{BI}=G(\boldsymbol{\sigma}_{BI})\boldsymbol{\omega}_{BI}^{B}$
(1a)
$G(\boldsymbol{\sigma}_{BI})=\dfrac{1}{2}\left(\dfrac{1-\boldsymbol{\sigma}^{T}_{BI}\boldsymbol{\sigma}_{BI}}{2}\textbf{{I}}_{3}+\boldsymbol{\sigma}^{\times}_{BI}+\boldsymbol{\sigma}_{BI}\boldsymbol{\sigma}^{T}_{BI}\right)$
$\textbf{{J}}\dot{\boldsymbol{\omega}}_{BI}^{B}+\boldsymbol{\omega}_{BI}^{B}\times\textbf{{J}}\boldsymbol{\omega}_{BI}^{B}=\boldsymbol{\tau}_{c}^{B}+\boldsymbol{\tau}_{d}^{B}$
(1b)

where $\boldsymbol{\omega}_{BI}^{B}\in\mathbb{R}^{3}$ denotes the angular
velocity expressed in the body-fixed frame,
$\textbf{{J}}\in\mathbb{R}^{3\times 3}$ is the inertia matrix,
$\textbf{{I}}_{3}$ denotes the identity matrix, and $(\textbf{{x}})^{\times}$
is the $3\times 3$ skew-symmetric cross-product matrix associated with vector
$\textbf{{x}}\in\mathbb{R}^{3}$. $\boldsymbol{\tau}_{c}^{B}$ and
$\boldsymbol{\tau}_{d}^{B}$ represent the control torque and the disturbance,
respectively.

$\boldsymbol{\sigma}_{XY}$ denotes the orientation of X frame relative to Y
frame. $\omega^{X}_{YZ}$ is the angular velocity of Y frame relative to Z
frame expressed in X frame. Then, the relative attitude between two frames is
defined as

$\boldsymbol{\sigma}_{\scriptscriptstyle
XY}=\dfrac{\boldsymbol{\sigma}_{\scriptscriptstyle
YI}(\boldsymbol{\sigma}^{T}_{\scriptscriptstyle
XI}\boldsymbol{\sigma}_{\scriptscriptstyle
XI}-1)+\boldsymbol{\sigma}_{\scriptscriptstyle
XI}(1-\boldsymbol{\sigma}^{T}_{\scriptscriptstyle
YI}\boldsymbol{\sigma}_{\scriptscriptstyle
YI})-2\boldsymbol{\sigma}_{\scriptscriptstyle
YI}^{\times}\boldsymbol{\sigma}_{\scriptscriptstyle
XI}}{1+\boldsymbol{\sigma}^{T}_{\scriptscriptstyle
XI}\boldsymbol{\sigma}_{\scriptscriptstyle
XI}\boldsymbol{\sigma}^{T}_{\scriptscriptstyle
YI}\boldsymbol{\sigma}_{\scriptscriptstyle
YI}+2\boldsymbol{\sigma}^{T}_{\scriptscriptstyle
YI}\boldsymbol{\sigma}_{\scriptscriptstyle XI}}$

and the dynamics of $\boldsymbol{\sigma}_{\scriptscriptstyle XY}$ is given by

$\dot{\boldsymbol{\sigma}}_{XY}=G(\boldsymbol{\sigma}_{XY})\boldsymbol{\omega}_{XY}^{X}$
(2a)
$\textbf{{J}}\dot{\boldsymbol{\omega}}_{XY}^{X}+\boldsymbol{\omega}_{XI}^{X}\times\textbf{{J}}\boldsymbol{\omega}_{XI}^{X}-\textbf{{J}}(\boldsymbol{\omega}^{X}_{XY}\times\boldsymbol{\omega}^{X}_{YI}-\textbf{{
C}}^{X}_{Y}\dot{\boldsymbol{\omega}}_{YI}^{Y})=\boldsymbol{\tau}_{c}^{B}$ (2b)

where
$\boldsymbol{\omega}^{X}_{XY}=\boldsymbol{\omega}^{X}_{XI}-\boldsymbol{\omega}^{X}_{YI}$
and $\boldsymbol{\omega}^{X}_{YI}=\textbf{{
C}}^{X}_{Y}\boldsymbol{\omega}^{Y}_{YI}$. The rotation matrix in terms of the
MRPs from $Y$ frame to $X$ frame can be expressed as

$\textbf{{
C}}^{X}_{Y}=\textbf{{I}}_{3}+\dfrac{8(\boldsymbol{\sigma}^{\times}_{XY})^{2}-4(1-\boldsymbol{\sigma}^{T}_{XY}\boldsymbol{\sigma}_{XY})\boldsymbol{\sigma}^{\times}_{XY}}{(1+\boldsymbol{\sigma}^{T}_{XY}\boldsymbol{\sigma}_{XY})^{2}}$
(3)

The following properties will be frequently used in this paper:

$\boldsymbol{\sigma}^{T}_{XY}G(\boldsymbol{\sigma_{XY}})=\left(\dfrac{1+\boldsymbol{\sigma}^{T}_{XY}\boldsymbol{\sigma_{XY}}}{4}\right)\boldsymbol{\sigma_{XY}}^{T}$
(4)
$G(\boldsymbol{\sigma}_{XY})^{T}G(\boldsymbol{\sigma})=\left(\dfrac{1+\boldsymbol{\sigma}^{T}\boldsymbol{\sigma}}{4}\right)^{2}\textbf{{I}}_{3}$
(5)

According to the description of MPRs in [33], MRPs have geometric
singularities when $\phi=\pm 360^{\circ}$, and it is not unique because of the
shadow set, i.e.,
$\boldsymbol{\sigma}=\boldsymbol{\sigma}^{s},\boldsymbol{\sigma}^{s}=-\boldsymbol{\sigma}/\boldsymbol{\sigma}^{T}\boldsymbol{\sigma}$.
Recalling the definition of $\boldsymbol{\sigma}$, one knows that
$\|\boldsymbol{\sigma}\|\leq 1$ for all $|\phi|\leq 180^{\circ}$. Thus, the
spacecraft attitude can be globally parameterized with the shortest principal
rotation by switching the $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}^{s}$
at the unit sphere $\|\boldsymbol{\sigma}\|=1$. Consequently, we stipulate
that the magnitude of $\boldsymbol{\sigma}$ is bounded by 1, i.e.,
$\|\boldsymbol{\sigma}\|\leq 1$, which is suited to describe any
reorientation.

###  2.2 State and control constraints

The pointing constraint, the angular velocity constraint and the input
limitation are considered in this paper. For the pointing constraint, we
suppose the instantaneous angle $\vartheta$ between a body-fixed unit vector
$\textbf{{r}}^{B}_{c}$ (such as cameras) and a inertial constant unit vector
$\textbf{{r}}^{I}_{t}$ (observed target) should be maintained in a half-cone
angel $\vartheta_{m}$, i.e., $\theta\leq\vartheta_{m}$, which is equivalent to

$\mathcal{C}_{p}=\left\\{(\boldsymbol{\sigma}_{BI},\boldsymbol{\omega}_{BI}^{B}):\textbf{{r}}^{B}_{c}\cdot\textbf{{r}}^{B}_{t}\geq{\rm
cos}(\vartheta_{m}),\vartheta_{m}\in(0,\frac{\pi}{2})\right\\}$ (6)

where $\textbf{{r}}^{B}_{t}=\textbf{{ C}}^{B}_{I}\textbf{{r}}^{I}_{t}$ is the
expression of $\textbf{{r}}^{I}_{t}$ in $\mathscr{F}_{B}$.

In consideration of the payload requirements, the angular velocity constraint
is always exists. Then the constraint set is given by

$\mathcal{C}_{\omega}=\left\\{(\boldsymbol{\sigma}_{BI},\boldsymbol{\omega}_{BI}^{B}):\|\boldsymbol{\omega}_{BI}^{B}\|\leq\omega_{\rm
max},~{}\omega_{\rm max}>0\right\\}$ (7)

where $\omega_{max}\in\mathbb{R}^{+}$ is the maximum angular velocity
amplitude.

The angular momentum exchange devices such as reaction wheels and control
moment gyros are usually used as the spacecraft attitude control actuators.
These devices may be saturated when the command torque is large. For
simplicity, the actuator constraint is formulated as

$\mathcal{C}_{\tau}=\left\\{(\boldsymbol{\sigma}_{BI},\boldsymbol{\omega}_{BI}^{B}):\|\boldsymbol{\tau}_{c}^{B}\|\leq\tau_{\rm
max},~{}\tau_{\rm max}>0\right\\}$ (8)

where $\tau_{\rm max}\in\mathbb{R}^{+}$ is the maximum allowable control
torque.

Finally, the dynamic safety margin of the system is the intersection of the
aforementioned three subsets:

$\mathcal{C}=\mathcal{C}_{p}\cap\mathcal{C}_{\omega}\cap\mathcal{C}_{\tau}$
(9)

###  2.3 Problem statements

This paper aims to develop a ERG control scheme that drives the system states
$(\sigma_{BI},\boldsymbol{\omega}_{BI}^{B})$ to the desired equilibrium
$(\sigma_{DI},\textbf{0}_{3\times 1})$ while satisfying the constraints (6),
(7), and (8). The proposed ERG-based control structure (shown in Fig. 1)
consists of two cascaded control units. The primary controller is given by an
angular velocity observer-based output feedback controller, which is able to
pre-stabilize the unconstrained system to an auxiliary reference
$\sigma_{VI}$. The reference governor (navigation layer) unit is designed to
guarantee the constraint enforcement by manipulating the kinematics of
$\sigma_{VI}$. Clearly, the asymptotic convergence property of the closed-loop
control system will be achieved by the goal that auxiliary reference
$\sigma_{VI}$ asymptotically tends to $\sigma_{DI}$. The reference attitude
$\sigma_{DI}$ remains constant and is chosen within the admissible region,
i.e., equation (6) holds true when $\sigma_{BI}=\sigma_{DI}$.

Figure 1:  The architecture of the explicit reference governor based attitude
control scheme.

##  3 Primary controller design

This section proposes an angular velocity-free control law so as to stabilize
the attitude toward a constant reference $\boldsymbol{\overline{\sigma}}_{VI}$
($\boldsymbol{\dot{\sigma}}_{VI}=\textbf{0}_{3\times 1}$) when the constraints
and disturbance are neglected ($\boldsymbol{\tau}_{d}=\textbf{0}_{3\times
1}$). The time variation of $\sigma_{VI}$ will be addressed by the reference
management unit, which is detailed in the next section.

###  3.1 Angular velocity observer design

The angular velocity observer is constructed based on the I&I theory [30, 31].
Let $\mathscr{F}_{E}$ be the estimation frame of the $\mathscr{F}_{B}$. The
attitude and angular estimation errors in terms of MRPs are given by

$\boldsymbol{\sigma}_{\scriptscriptstyle
BE}=\dfrac{\boldsymbol{\sigma}_{\scriptscriptstyle
EI}(\boldsymbol{\sigma}^{T}_{\scriptscriptstyle
BI}\boldsymbol{\sigma}_{\scriptscriptstyle
BI}-1)+\boldsymbol{\sigma}_{\scriptscriptstyle
BI}(1-\boldsymbol{\sigma}^{T}_{\scriptscriptstyle
EI}\boldsymbol{\sigma}_{\scriptscriptstyle
EI})-2\boldsymbol{\sigma}_{\scriptscriptstyle
EI}^{\times}\boldsymbol{\sigma}_{\scriptscriptstyle
BI}}{1+\boldsymbol{\sigma}^{T}_{\scriptscriptstyle
BI}\boldsymbol{\sigma}_{\scriptscriptstyle
BI}\boldsymbol{\sigma}^{T}_{\scriptscriptstyle
EI}\boldsymbol{\sigma}_{\scriptscriptstyle
EI}+2\boldsymbol{\sigma}^{T}_{\scriptscriptstyle
EI}\boldsymbol{\sigma}_{\scriptscriptstyle BI}}$ (10a)
$\boldsymbol{\omega}_{BE}^{B}=\boldsymbol{\omega}_{BI}^{B}-C^{B}_{E}\boldsymbol{\omega}_{EI}^{E}=\boldsymbol{\omega}_{BI}^{B}-\boldsymbol{\omega}_{EI}^{B}$
(10b)

To ensure $\boldsymbol{\omega}_{BE}^{B}\rightarrow 0$ and
$\boldsymbol{\sigma}_{BE}\rightarrow 0$, a scalar $\varpi\in\mathbb{R}^{+}$
that can ’cover’ the $\boldsymbol{\omega}_{EI}^{B}$ is introduced as

$\varpi=\sqrt{\varepsilon_{\omega}+\|\boldsymbol{\omega}_{EI}^{B}\|^{2}}$ (11)

where $\varepsilon_{\omega}\in\mathbb{R}^{+}$ is a constant to be selected,
which is utilized to ensure the existence of the time derivative of $\varpi$.
Then, $\boldsymbol{\omega}_{EI}^{B}$ is generated by

$\boldsymbol{\omega}_{EI}^{B}=\boldsymbol{\xi}+4\textbf{{J}}^{-1}\beta(\underline{\varpi})\boldsymbol{\sigma}_{BE}$
(12)

where $\underline{\varpi}$ is the estimate of $\varpi$, $\boldsymbol{\xi}$ and
$\beta(\underline{\varpi})$ are the parameter related to
$\boldsymbol{\omega}_{EI}^{B}$ and a function of $\underline{\varpi}$,
respectively. The dynamics of $\boldsymbol{\xi}$, $\underline{\varpi}$, and
$\boldsymbol{\sigma}_{EI}$ are designed as

$\displaystyle\dot{\boldsymbol{\xi}}=$
$\displaystyle\textbf{{J}}^{-1}(-\boldsymbol{\omega}_{EI}^{B}\times\textbf{{J}}\boldsymbol{\omega}_{EI}^{B}+\boldsymbol{\tau}_{c}^{B})-4\textbf{{J}}^{-1}\dot{\beta}(\underline{\varpi})\boldsymbol{\sigma}_{BE}$
$\displaystyle-4\textbf{{J}}^{-1}\beta(\underline{\varpi})\dot{\boldsymbol{\sigma}}_{1BE}$
(13a) $\displaystyle\dot{\underline{\varpi}}=\varpi^{-1}$
$\displaystyle(\boldsymbol{\omega}_{EI}^{B})^{T}\textbf{{J}}^{-1}(-\boldsymbol{\omega}_{EI}^{B}\times\textbf{{J}}\boldsymbol{\omega}_{EI}^{B}+\boldsymbol{\tau}_{c}^{B})-K_{\underline{\varpi}}(\underline{\varpi}-\varpi)$
(13b)
$\displaystyle\dot{\boldsymbol{\sigma}}_{EI}=G(\boldsymbol{\sigma}_{EI})(\boldsymbol{\omega}_{EI}^{E}+K_{\sigma}\textbf{{C}}^{E}_{B}\boldsymbol{\sigma}_{BE})$
(13c)

where $K_{\underline{\varpi}}$ and $K_{\sigma}$ are the dynamic gains to be
designed, $\dot{\boldsymbol{\sigma}}_{1BE}$ represents part of the dynamics of
$\boldsymbol{\sigma}_{BE}$, and it can be obtained from (1a), (10b), and (13c)
:

$\displaystyle\dot{\boldsymbol{\sigma}}_{BE}=$ $\displaystyle
G(\boldsymbol{\sigma}_{BE})(\boldsymbol{\omega}_{BI}^{B}-\boldsymbol{\omega}_{EI}^{B})$
$\displaystyle=$ $\displaystyle
G(\boldsymbol{\sigma}_{BE})[\boldsymbol{\omega}_{BI}^{B}-\textbf{{C}}^{B}_{E}(\boldsymbol{\omega}_{EI}^{E}+K_{\sigma}\textbf{{C}}^{E}_{B}\boldsymbol{\sigma}_{BE})]$
$\displaystyle=$ $\displaystyle
G(\boldsymbol{\sigma}_{BE})(\boldsymbol{\omega}_{BE}^{B}-K_{\sigma}\boldsymbol{\sigma}_{BE})$
$\displaystyle=$
$\displaystyle\dot{\boldsymbol{\sigma}}_{1BE}+\dot{\boldsymbol{\sigma}}_{2BE}$
(14)

with
$\dot{\boldsymbol{\sigma}}_{1BE}=-G(\boldsymbol{\sigma}_{BE})K_{\sigma}\boldsymbol{\sigma}_{BE}$
and
$\dot{\boldsymbol{\sigma}}_{2BE}=G(\boldsymbol{\sigma}_{BE})\boldsymbol{\omega}_{BE}^{B}$.

To inject the nonlinear terms in the dynamics of
$\boldsymbol{\omega}_{BE}^{B}$ a dynamic scaling technique is introduced:

$\textbf{{z}}=\dfrac{\boldsymbol{\omega}_{BE}^{B}}{\textit{r}}$ (15)

where $r$ is the dynamic scaling factor and is updated by the following law

$\dot{\textit{r}}=\dfrac{\textit{r}}{J_{m}}(J_{M}\|\varpi-\underline{\varpi}\|)-\dfrac{k_{r}}{J_{M}}(r-1)$
(16)

where $k_{r}\in\mathbb{R}^{+}$ is the dynamic scaling gain to be determined,
$J_{m}$ and $J_{M}$ are the minimum and the maximum eigenvalues of the inertia
matrix J, respectively. If $r(t)=1$, $\dot{\textit{r}}\geq 0$. Hence, it
satisfies $r(t)\geq 1$ for all $t\geq 0$ when $r(0)\geq 1$. Finally, the
convergence analysis of the proposed observer (12) is summarized as the
following proposition.

######  proposition 1.

Consider the angular velocity observer in (12) with dynamics given in (13),
(16), and the gains are given as

$\beta(\underline{\varpi})=4\underline{\beta}(\underline{\varpi})G^{T}(\boldsymbol{\sigma}_{BE})$
(17a)
$\underline{\beta}(\underline{\varpi})=J_{M}\|\underline{\varpi}\|+\dfrac{J_{m}k_{r}}{J_{M}}+1+\rho_{\varpi}$
(17b)
$K_{\underline{\varpi}}=8\left(\dfrac{\|\boldsymbol{\omega}_{EI}^{B}\|\beta(\underline{\varpi})r}{J_{m}}\right)^{2}+\frac{1}{2}r^{2}J_{M}+\rho_{\underline{\varpi}}$
(17c) $K_{\sigma}=\frac{1}{2}r^{2}+\rho_{\sigma}$ (17d)
$k_{r}=\frac{1}{2}\dfrac{J^{2}_{M}}{J_{m}}+\rho_{r}$ (17e)

where $\rho_{\varpi}$, $\rho_{\underline{\varpi}}$, $\rho_{\sigma}$, and
$\rho_{r}$ are positive constants that can be tuned for different convergence
rates of the estimation errors. Then, the dynamic scaling factor $r$ is
bounded and the errors globally exponentially converges to the origin, i.e.,
${\rm lim}_{t\rightarrow\infty}e^{\alpha
t}\|\boldsymbol{\omega}_{BE}^{B}\|=0,\alpha\in\mathbb{R}^{+}$.

Proof: See the Appendix. $\hfill\blacksquare$

###  3.2 Velocity-free controller design

The aforementioned angular velocity observer is used to derive a velocity-free
feedback attitude controller. As shown in Fig. 1, the following theorem
summarize the result on the unconstrained output controller.

######  Theorem 1.

Consider the attitude dynamics given in (1) and the angular velocity observer
given in (12)-(17). Then, the output feedback control law is given by

$\boldsymbol{\tau}_{c}^{B}=-k_{p}\boldsymbol{\sigma}_{BV}-k_{d}\boldsymbol{\omega}_{EI}^{B}$
(18)

with $k_{p},~{}k_{d}>0$, the equilibrium
$(\boldsymbol{\overline{\sigma}}_{VI},\textbf{0}_{3\times 1})$ is
asymptotically stable within the admissible set, i.e., ${\rm
lim}_{t\rightarrow\infty}(\boldsymbol{\sigma}_{BI},\boldsymbol{\omega}_{BI}^{B})=(\boldsymbol{\overline{\sigma}}_{VI},\textbf{0}_{3\times
1})$.

Proof: By employing (10b), the control law (18) can be expressed as a full-
state controller augmented by perturbations resulting from velocity estimation
errors, namely,

$\boldsymbol{\tau}_{c}^{B}=-k_{p}\boldsymbol{\sigma}_{BV}-k_{d}\boldsymbol{\omega}_{BI}^{B}+k_{d}\boldsymbol{\omega}_{BE}^{B}$
(19)

Consider a Lyapunov function candidate as follows:

$V_{c}=2k_{p}{\rm
ln}(1+\boldsymbol{\sigma}^{2}_{BV})+\frac{1}{2}(\boldsymbol{\omega}_{BI}^{B})^{T}\textbf{{J}}\boldsymbol{\omega}_{BI}^{B}$
(20)

Taking the time derivative of (20) along (1), (4), and (19), one can obtain

$\displaystyle\dot{V}_{c}=$ $\displaystyle
4k_{p}\dfrac{\boldsymbol{\sigma}_{BV}^{T}\boldsymbol{\dot{\sigma}}_{BV}}{1+\boldsymbol{\sigma}_{BV}^{2}}+(\boldsymbol{\omega}_{BI}^{B})^{T}\textbf{{J}}\boldsymbol{\dot{\omega}}_{BI}^{B}$
$\displaystyle\leq$ $\displaystyle-
k_{d}\|\boldsymbol{\omega}_{BI}^{B}\|^{2}+k_{d}\|\boldsymbol{\omega}_{BI}^{B}\|\|\boldsymbol{\omega}_{BE}^{B}\|$
(21)

Clearly, $\dot{V}_{c}$ contains a term with indefinite sign induced by angular
velocity estimation error. To eliminate this effect, consider a positive
definite Lyapunov function in the following form:

$V=V_{c}+\delta_{z}V_{z}$ (22)

where $\delta_{z}$ is a positive constant to be determined. Differentiating
$V$ and applying (21) and (39) yields

$\displaystyle\dot{V}\leq$ $\displaystyle-
k_{d}\|\boldsymbol{\omega}_{BI}^{B}\|^{2}+k_{d}\|\boldsymbol{\omega}_{BI}^{B}\|\|\boldsymbol{\omega}_{BE}^{B}\|-\delta_{z}(1+\rho_{\varpi})\|\textbf{{z}}\|^{2}$
(23) $\displaystyle\leq$
$\displaystyle-[\|\boldsymbol{\omega}_{BI}^{B}\|\|\boldsymbol{\omega}_{BE}^{B}\|]\\!\\!\begin{bmatrix}k_{d}\\!\\!\\!&-\dfrac{1}{2}k_{d}\\\
-\dfrac{1}{2}k_{d}\\!\\!\\!&\delta_{z}r^{-2}(1+\rho_{\varpi})\end{bmatrix}\\!\\!\begin{bmatrix}\|\boldsymbol{\omega}_{BI}^{B}\|\\\
\|\boldsymbol{\omega}_{BE}^{B}\|\end{bmatrix}\\!$

Given that $1\leq r<\infty$, there exists a sufficiently large $\delta_{z}$
such that $\dot{V}$ is negative semi-definite for
$\boldsymbol{\overline{\sigma}}{VI}$. By using the LaSalle invariance
principle, one can conclude that the equilibrium point
$(\boldsymbol{\overline{\sigma}}_{VI},\textbf{0}_{3\times 1})$ of the system
is asymptotically stable. This completes the proof. $\hfill\blacksquare$

Obviously, $\boldsymbol{\sigma}_{VI}$ is time varying, Theorem 1 addresses the
claim concerning the tracking error stability of
$\boldsymbol{\overline{\sigma}}_{VI}$ rather than $\boldsymbol{\sigma}_{VI}$.
In fact, since the final state $\boldsymbol{\sigma}_{DI}$ is a constant
attitude, the inner loop controller only needs to ensure that the attitude can
converge to the final state. Additionally, the precise acquisition of angular
velocity is challenging, resulting in difficulties in strictly guaranteeing
the angular velocity constraint. Fortunately, the value of the dynamic scaling
factor $r$ reflects the estimation error. By designing functions associated
with $r$, it becomes possible to satisfy the angular velocity constraints.
These properties will be utilized in the subsequent section.

######  Remark 1.

Throughout the preceding analysis, it is evident that the design of the output
feedback controller remains independent of the angular velocity observer
(refer to (18) and (19)), which greatly reduces the difficulty of the
controller design. Moreover, (13) and (17) suggest a conflict between the
scale of the observer gains and system robustness. However, there are always
parameter uncertainties in practical missions, so it is necessary to balance
the two properties, which will be verified in detail in the simulation.

##  4 Reference management

The reference management layer of ERG (illustrated in Fig. 1) devises an
auxiliary control law that manipulates the reference state towards the primary
stabilized system [36, 20]. The objective of this section is to handle the
constraints outlined in (6) $-$ (8) by designing the safety margin and the
navigation field, which is achieved by the invariant set in the Lyapunov
function centered on the reference state $\sigma_{VI}$. The formulation of the
auxiliary reference takes the following form:

$\dot{\sigma}_{VI}=\Delta(\sigma_{BV},\boldsymbol{\omega}_{EI}^{B})\chi(\sigma_{VI},\sigma_{VD})$
(24)

where
$\Delta(\sigma_{BV},\boldsymbol{\omega}_{EI}^{B}):\mathbb{R}^{3}\times\mathbb{R}^{3}\rightarrow\mathbb{R}$
is the dynamic safety margin that indicates how safe it is within the
allowable set.
$\chi(\sigma_{VI},\sigma_{VD}):\mathbb{R}^{3}\times\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$
denotes the navigation field of the current state $\sigma_{VI}$, and the
$\sigma_{VD}$ is utilized to drive the $\sigma_{VI}$ towards to $\sigma_{DI}$.

###  4.1 Safety margin

From an intuitive perspective, the safety margin can be treated as the
distance between the constraint boundary and the navigation field. Since
$\dot{V}$ is negative semi-define (see (23)), the forward invariant set
$\left\\{(\boldsymbol{\sigma}_{BI},\boldsymbol{\omega}_{BI}^{B}):V\leq\Gamma\right\\}$
can be used to design the safety margin, where the upper bound
$\Gamma(\boldsymbol{\sigma}_{VI}^{V},\boldsymbol{\omega}_{VI}^{V})$ is
determined by the constraints (6) $-$ (8). In [21] and [22], authors design
the dynamic safety margin in the form
$\Delta(\sigma_{BV},\boldsymbol{\omega}_{EI}^{B})=k_{e}(\Gamma-V)$, where the
constant $k_{e}$ is used to adjust the dynamic performance. Unfortunately,
since the exact estimation error $\boldsymbol{\omega}_{BE}^{B}$ is
unavailable, the angular velocity can not be obtained either. In order to
prevent $\Delta$ being negative caused by $\boldsymbol{\omega}_{BE}^{B}$,
$\Delta$ can be designed as

$\Delta(\sigma_{BV},\boldsymbol{\omega}_{EI}^{B})=\left\\{\begin{aligned}
k_{e}(\Gamma-V),~{}&~{}\Gamma>V\\\ 0~{}~{}~{}~{}~{}~{}~{},~{}&~{}\Gamma\leq
V\end{aligned}\right.$ (25)

####  4.1.1 Pointing constraint

The geometric relationship about the pointing constraint is displayed in Fig.
2, where $\vartheta_{e}\in(0,\frac{\pi}{2})$ is the safety margin of
$\sigma_{BV}^{B}$ and satisfies $\vartheta_{e}=\vartheta_{m}-\vartheta$. When
the body frame is coincided with the reference frame $\mathscr{F}_{V}$, the
pointing angle $\overline{\vartheta}$ is denote by

$\overline{\vartheta}=arcos(\overline{\textbf{{ r}}}^{V}_{c}\cdot\textbf{{
r}}^{V}_{t})$ (26)

Note that $\overline{\textbf{{ r}}}^{V}_{c}=\textbf{{ r}}^{B}_{c}$ is a
virtual constant unit vector expressed in $\mathscr{F}_{V}$ rather than
$\overline{\textbf{{ r}}}^{V}_{c}=\textbf{{C}}^{V}_{B}\textbf{{ r}}^{B}_{c}$.
Under these conditions, the safety margin of $\sigma_{VI}^{V}$ satisfies
$\overline{\vartheta}_{e}=\vartheta_{m}-\overline{\vartheta}$ and
$\overline{\vartheta}_{e}\in(0,\frac{\pi}{2}]$. Let $\alpha\in[0,\pi]$ be the
gap between $\textit{{n}}_{BV}$ and the unit vector $\textbf{{ r}}^{B}_{c}$,
then it satisfies

$\textbf{{ r}}^{B}_{c}\cdot\textit{{n}}_{BV}={\rm cos}(\alpha)$ (27)

Figure 2: Constrained attitude region.

Let $\varphi\in(0,\pi)$ denote the orientation from $\textbf{{ r}}^{B}_{c}$ to
$\overline{\textbf{{ r}}}^{V}_{c}$,and they have the following relationship:

${\rm sin}(\dfrac{\varphi}{2})={\rm sin}(\dfrac{\phi_{BV}}{2}){\rm
sin}(\alpha)$ (28)

Obviously, there is a positive correlation between $\varphi$ and $\phi_{BV}$,
and they satisfies $\varphi\leq\phi_{BV}$. Since
$|\vartheta_{e}-\overline{\vartheta}_{e}|\leq\varphi$, if
$\varphi\leq\overline{\vartheta}_{e}$, then $\vartheta_{e}\geq 0$ can be
guaranteed.

According to [7], if $\boldsymbol{\omega}_{EI}^{E}=\textbf{0}_{3\times 1}$
($\boldsymbol{\omega}_{BI}^{B}$ is precisely known), then
$\dot{V}_{c}=\dot{V}\leq 0$, the threshold of $V\leq\Gamma_{p}^{\prime}$ can
be designed as

$\Gamma_{p}^{\prime}=\left\\{\begin{aligned} 2k_{p}{\rm
ln}\left\\{1+\left(\frac{1-\sqrt{1-a_{p}^{2}}}{a_{p}}\right)^{2}\right\\}&~{}~{}\alpha\in(0,\pi)\\\
\infty~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}&~{}~{}\alpha=0~{}{\rm
or}~{}\pi\end{aligned}\right.$ (29)

where $a_{p}={\rm sin}(\dfrac{\overline{\vartheta}_{e}}{2})/{\rm
sin}(\alpha)$.

When $V=\Gamma_{p}^{\prime}$, $\dot{\sigma}_{VI}=\textbf{0}_{3\times 1}$.
According to Theorem1, the time derivate of $V$ and $\Gamma_{p}^{\prime}$
satisfies $\dot{V}\leq 0,{\dot{\Gamma}}^{\prime}_{p}=0$, which means the
pointing constraint (6) will never be violated. However,
$\boldsymbol{\omega}_{BI}^{B}$ and $\delta_{z}$ are unavailable, which means
$V$ and $V_{c}$ are unavailable, thus they can not be used to design the
threshold of the pointing constraint. Hence, we approximate the Lyapunov
function V with the following form

$\overline{V}=2k_{p}{\rm
ln}(1+\boldsymbol{\sigma}_{BV})+\frac{1}{2}(\boldsymbol{\omega}_{EI}^{B})^{T}\textbf{{J}}\boldsymbol{\omega}_{EI}^{B}$
(30)

Although $\dot{\overline{V}}$ is sign indefinite, according Theorem 1,
$\overline{V}$ is asymptotically convergent, i.e., ${\rm
lim}_{t\rightarrow\infty}\overline{V}(t)=0$. Accordingly, the threshold of the
pointing constraint is designed as

$\Gamma_{p}=\dfrac{\Gamma_{p}^{\prime}}{r^{k_{1}}}$ (31)

where $k_{1}>0$ is a constant parameter. Similarly to the analysis in Sec.
3.1, the larger $\boldsymbol{\omega}_{BE}^{B}$ is, the lager $r$ is, and the
smaller $\Gamma_{p}$ is. Although the exact relationship between
$\boldsymbol{\omega}_{BE}^{B}$ and $r$ is unknown, by tuning $k_{1}$, a
conservative but safe threshold of the pointing constraint without angular
velocity measurement can be obtained.

####  4.1.2 Angular velocity constraint

The angular velocity constraint given in (7) is a convex set. As discussed in
Sec. 4.1.1, when $\boldsymbol{\omega}_{EI}^{E}=\textbf{0}_{3\times 1}$, the
threshold of the angular velocity constraint $\Gamma_{\omega}$ can be selected
as

$\Gamma^{\prime}_{\omega}=\frac{1}{2}J_{m}\omega_{\rm max}^{2}$ (32)

Similar to the pointing constraint, when $V=\Gamma^{\prime}_{\omega}$,
$\dot{\sigma}_{VI}=\textbf{0}_{3\times 1}$. Recalling Theorem1, the time
derivate of them satisfies $\dot{V}\leq 0,{\dot{\Gamma}}^{\prime}_{\omega}=0$,
which means the angular velocity constraint (7) will never be violated. When
$\boldsymbol{\omega}_{EI}^{E}\neq\textbf{0}_{3\times 1}$, the Lyapunov
function is approximate by (30), and the threshold can be selected as

$\Gamma_{\omega}=\frac{\Gamma^{\prime}_{\omega}}{r^{k_{2}}}$ (33)

where $k_{2}>0$ is a constant parameter used to tuning $\Gamma_{\omega}$.

####  4.1.3 Actuator saturation

Similar to the preceding constraints, $\boldsymbol{\omega}_{BI}^{B}$ is
substituted with $\boldsymbol{\omega}_{EI}^{B}$ and we omit the estimation
error. Adhering to the approach outlined in [22], the saturation constraint
(8) can be satisfied by solving the subsequent optimization problem

Problem

${\rm min}~{}2k_{p}{\rm
ln}(1+\boldsymbol{\sigma}_{BV})+\frac{1}{2}(\boldsymbol{\omega}_{EI}^{B})^{T}\textbf{{J}}\boldsymbol{\omega}_{EI}^{B}$

subject to

$\lvert\boldsymbol{\sigma}_{BV}\rvert_{i}\leq 1$ (34a) $\lvert
k_{p}\boldsymbol{\sigma}_{BV}+k_{d}\boldsymbol{\omega}_{EI}^{B}\rvert_{i}\geq\tau_{\rm
max}$ (34b)

Then the threshold $\Gamma_{\tau}$ can be obtained by selecting the minimum
value from the aforementioned optimization problem for $i={1,2,3}$.
Consequently, the upper-bound of the system subject to the constraint (9) can
be concluded as $\Gamma={\rm
min}\\{\Gamma_{p},\Gamma_{\omega},\Gamma_{\tau}\\}$, a conclusion that can be
substantiated by employing the same arguments presented in [21].

###  4.2 Navigation layer

The navigation field $\chi(\sigma_{VI},\sigma_{VD})$ will be designed in this
section to ensure that the auxiliary reference $\boldsymbol{\sigma}_{BV}$
towards to the desired reference $\boldsymbol{\sigma}_{BD}$. Consequently, the
trajectory of $\chi(\sigma_{VD})$ must strictly align within the permissible
set $\mathcal{C}$. Given that both the initial and final attitudes fall within
the constraints and the pointing constraint $\mathcal{C}_{p}$ constitutes a
convex set, the shortest distance on the attitude manifold adheres to the
constraints. The navigation trajectory $\chi(\sigma_{VD})$ is designed by

$\chi(\sigma_{VD})=-G(\boldsymbol{\sigma}_{VD})\boldsymbol{\sigma}_{VD}$ (35)

Since $\boldsymbol{\omega}_{BI}^{B}=\textbf{0}_{3\times 1}$ and
$\boldsymbol{\tau}_{c}^{B}=\textbf{0}_{3\times 1}$ represent the equilibrium
point, the constraints (7) and (8) are always satisfied at steady-state. Then,
the main results about the constrained attitude maneuver control without
angular velocity measurement is presented in the following proposition.

######  proposition 2.

Given the spacecraft attitude dynamics (1) subject to the constraints (9) with
the angular velocity observer (12) controlled by (18), and let (24) be the
navigation layer subject to the dynamic safety margin (25), and the navigation
field (35). Then, for any initial states satisfy the constraints and
$V(0)\leq\Gamma(0)$, the following statements hold. 1) For any constant
reference $\boldsymbol{\sigma}_{DI}\in\mathcal{C}$, the system constraints are
all satisfied. 2) The auxiliary reference $\boldsymbol{\sigma}_{VI}$ updated
by (24) asymptotically converges to $\boldsymbol{\sigma}_{DI}$.

Proof: See [7]. $\hfill\blacksquare$

######  Remark 2.

By integrating the ERG design process with the aforementioned analysis, it can
be seen that by designing a trajectory from the current state to the final
state that satisfies the constraint conditions, and then the controller drives
the system state and the reference state error within a certain range, the
system states can be guaranteed to reach the target state while the
constraints are met. An additional advantage of this strategy lies in its
capability to maintain effective control even in the absence of state
constraints. In comparison to alternative control algorithms like PID, the ERG
algorithm can track the ${\sigma}_{VD}$ independently generated by the
reference management with smaller error than ${\sigma}_{BD}$. This capability
leads to improved control performance characterized by enhanced speed and
precision.

##  5 Numerical simulations

This section aims to demonstrate the effectiveness of the proposed angular
velocity free attitude control algorithm in the presence of multi-constraints.
The objective involves maneuvering the rigid spacecraft from a specific
initial state to a pre-defined target, incorporating considerations for
attitude constraints, angular velocity constraints, and control saturation
concurrently. Besides, Monte Carlo results are conducted to further verify the
robustness of the proposed control scheme. The inertia of the spacecraft is
given by

$\textbf{\emph{J}}=\begin{bmatrix}15.2&-1&2\\\ -1&18.3&-0.5\\\
2&-0.5&16.1\end{bmatrix}\mathrm{kg.m}^{2}.$

The initial states are set as
$\boldsymbol{\sigma}_{BI}(0)=[-0.119,0.000,0.159]^{T}$ and
$\boldsymbol{\omega}(0)=[0,-0.01,0.01]^{T}{\rm rad/s}$. The constraint
conditions and target state are chosen in Table 1. Besides, the threshold of
actuator saturation $\Gamma_{\tau}$ is obtained by solving from the $Problem$
via $fmincon$ function in Matlab 2021, which is 0.0468, and the observer
parameters are shown in Table 2. The control elements are selected as
$k_{p}=1.5$ and $k_{d}=2.5$. For the brevity and intuition, Euler angles
$[\varphi,\vartheta,\psi]^{T}$ with sequence $3-1-2$ are used to plot the
attitude.

Table 1: System Constraint Conditions. Parameters | Values
---|---
$\boldsymbol{\sigma}_{DI}$ | $[0,0,0]^{T}$
$\vartheta_{m}$ | $38^{\circ}$
$\textbf{{r}}^{I}_{t}$ | $[1/\sqrt{3},-1/\sqrt{3},1/\sqrt{3}]^{T}$
$\textbf{{r}}^{B}_{c}$ | $[0,-1/\sqrt{2},1/\sqrt{2}]^{T}$
$\omega_{max}$ | $0.035{\rm rad/s}$
$\tau_{\rm max}$ |  $0.1$ N.m.s
$k_{e}$ | $1000$
$k_{1},k_{2}$ | $2$
Table 2: Observer parameters. Parameters | Values
---|---
$J_{m}$ | $18.3$
$J_{M}$ | $15.2$
$\rho_{\sigma},\rho_{\varpi},\rho_{\underline{\varpi}},\rho_{r},\varepsilon_{\omega}$ | $0.1$
$r(0)$ | $1$
$\boldsymbol{\xi}(0)$ | $[0,0,0]^{T}$
$\boldsymbol{\sigma}_{EI}(0)$ | $[-0.119,0.000,0.159]^{T}$
$\boldsymbol{\omega}_{EI}^{B}$ | $[0,0,0]^{T}{\rm rad/s}$

###  5.1 Performance of the proposed control scheme

The simulation results are depicted in Figs. 4 $-$ 12, where the dash curves
illustrate the simulations conducted without reference management i.e., the
constraints are not integrated into the control scheme. The attitude
trajectories and angular velocities depicted in Figs. 4 and 4 indicate that
with the navigation layer, the trajectories of attitude and the velocity
become smoother and the overshoot is smaller.

Figure 3: Attitude trajectories.

Figure 4: Angular velocities.

Figure 5:  Attitude and velocity estimation errors.

Figure 6: Injection gain r.

Figure 7: Pointing constraint.

Figure 8: Angular velocity constraint.

Fig. 6 illustrates the attitude and angular velocity estimation errors
generated by the I$\&$I based observer designed in (12)$-$(17). In the
logarithmic scale, the estimation errors $\|\boldsymbol{\omega}_{BE}^{B}\|$
and $\|\boldsymbol{\sigma}_{BE}\|$ decreases in an almost straight line,
indicating that the estimation errors are exponentially convergent. An
intriguing observation arises: unlike $\|\boldsymbol{\omega}_{BE}^{B}\|$
decreases consistently over time, $\|\boldsymbol{\sigma}_{BE}\|$ initially
remains very small but undergoes an initial increase before subsequently
decreasing. This behavior arises because $\|\boldsymbol{\sigma}_{BE}\|$ is
utilized as an ”indicator” to assess the adequacy of the estimation of
$\|\boldsymbol{\omega}_{BE}^{B}\|$. As we set $\|\boldsymbol{\sigma}_{BE}\|=0$
as the initial condition, and the estimation error
$\|\boldsymbol{\omega}_{BE}^{B}(0)\|$ is substantial,
$\|\boldsymbol{\sigma}_{BE}\|$ grows initially. As
$\|\boldsymbol{\omega}_{BE}^{B}\|$ diminishes, $\|\boldsymbol{\sigma}_{BE}\|$
subsequently adjusts accordingly. As we can also seen from Figs. 6 and 6, due
to the large estimation error at the beginning, the injection gain $r$ is also
relatively large, but as the estimation error diminishes, $r$ swiftly
converges towards 1. These outcomes underscore the pivotal role of the
injection gain $r$ in the observer system, contributing significantly to
achieving the desired effectiveness of the observer. Furthermore, the final
angular velocity estimation error $\|\boldsymbol{\omega}_{BE}^{B}(0)\|$
measures approximately $10^{-4}$, which is mainly restricted by the simulation
setting 0.01s. If the step size is further reduced,
$\|\boldsymbol{\omega}_{BE}^{B}(0)\|$ can also be reduced.

Figure 9: Control torque limitation.

Figure 10: Reference trajectories.

Figure 11: Threshold value.

Figure 12: Attitude trajectories.

The pointing constraint, angular velocity constraint and the control torque
limitation are plotted in Figs. 8 $-$ 10. Evidently, in the absence of
reference management, the pointing angle exceeds the boundary at 23 seconds,
the angular velocity surpasses the limitation around 15 seconds, and the
actuator saturation occurs within the initial 10 seconds. Upon the
implementation of the navigation layer within the system, all three
constraints remain within permissible bounds. This is because the tracking
error can always maintains a small error relative to the reference
trajectories (see Fig. 10). Analyzing Figs. 10 and 12, reveals that the
convergence speed of the reference trajectory aligns closely with the
threshold error $\Gamma-V$. Overall, the simulation results are in line with
the theoretical analysis, verifying the effectiveness of the algorithm.

###  5.2 Monte Carlo simulation under disturbances

The aforementioned section presented numerical simulations without any
disturbance i.e., $\boldsymbol{\tau}_{d}=\textbf{0}_{3\times 1}$. This section
conducts Monte Carlo simulations encompassing disturbances to showcase the
robustness of the proposed ERG and I&I technology-based attitude control
scheme. To conduct the Monte Carlo simulations, randomized initial conditions
and parameters are detailed in Table 3. These newly introduced initial states
and parameters are integrated with previously selected simulation conditions
mentioned in the previous section, amounting to a total of 200 Monte Carlo
simulations. In order to ensure that all the cases can reach the final states,
those cases that do not meet the constraints at the initial moment will be
excluded, and the simulations last for 150 seconds. Besides, the external
disturbances are given as follows

$\boldsymbol{\tau}_{d}^{I}=\begin{pmatrix}\begin{bmatrix}2\\\ -1\\\
-3\end{bmatrix}+\begin{bmatrix}0.4\mathrm{sin}(\omega_{A}t+1.6)\\\
2\mathrm{sin}(\omega_{A}t+1.1)\\\
0.7\mathrm{sin}(\omega_{A}t-2.1)\end{bmatrix}\end{pmatrix}\times
10^{-5}\mathrm{N.m}^{2}$

where $\omega_{At}=0.01rad/s$.

Figs. 12 and 14 depict the attitude and angular velocity trajectories of Monte
Carlo simulations. Despite external disturbances and angular velocity
estimation errors influencing the control performance, the spacecraft achieved
the desired state smoothly. It can be seen from (12) and (13) that the
precision of the input torque significantly influences estimation errors. In
the Monte Carlo simulations, the presence of unknown disturbances results in
slightly larger estimation errors, as depicted in Figs. 14 and 16, compared to
those in Fig. 6. Nevertheless, these errors remain comparatively small,
affirming the observer’s robustness as designed in (12). Besides, the three
constraints are basically satisfied during the simulation (see Figs. 16 $-$
18). Meanwhile, one can also see that the angular velocity observer is
independent of the controller. No matter what constraints the system needs to
meet, the observer can converge without being affected by them. This
intriguing property enables us to enhance the ERG-based constrained controller
without limitations imposed by the observer.

Table 3: Randomized initial states and parameters. Variables | Ranges
---|---
$\textbf{{n}}(0)$ | $[-0.7,-0.5]\times[-0.1,0.1]\times[0.7,0.9]$
$\phi(0)$, rad | $[0.15\pi,0.35\pi]$
$\boldsymbol{\omega}(0)$, rad/s | $\\{[-2,2]\times[-2,2]\times[-2,2]\\}\times 10^{-3}$
$J_{m}$ | [13, 17]
$J_{M}$ | [15, 21]
$k_{p}=1.5$ | [1, 1.5]
$k_{d}=2.5$ | [2.5, 3.0]
$k_{e}=1000$ | [900, 1100]

Figure 13: Angular velocity trajectories.

Figure 14: Angular velocity estimation errors.

However, a slight setback is identified as there are two cases of slow
convergence (as denoted by the red circles in Figs. 12, 14, 16, and 18). This
occurrence is due to these specific cases being at the verge of the pointing
constraint initially. Due to external disturbances and angular velocity
estimation errors, the pointing angle slightly overflows the boundary.
Fortunately, the algorithm designed in this study addresses these effects by
employing an enhanced safety margin as described in (25), swiftly reeling back
the pointing angle. To prevent such occurrences rigorously, adjusting the
reference state’s margin suffices, a topic to be extensively explored in our
subsequent research. In spite of this, the proposed ERG-based constrained
controller still accomplished the maneuver objective with remarkable accuracy,
which in turn demonstrates its robustness against uncertainties.

Figure 15: Attitude estimation errors.

Figure 16: Distributions of pointing angles.

Figure 17: Distributions of angular velocity constraints.

Figure 18: Distributions of control torque.

##  6 Conclusion

This paper develop a constrained output feedback attitude reorientation
problem via ERG and I$\&$I technologies, where pointing constraint, angular
velocity constraint, and the control saturation are considered. The stability
of a angular velocity observer and the output feedback controller is roughly
proved. The inner loop of ERG is conducted by the angular velocity observer-
based output feedback controller, and the navigation layer is designed by
manipulating the auxiliary reference state without violating the constraints
while asymptotically converges to the desired reference. The performance of
the proposed angular velocity observer and the ERG is meticulously analyzed
and discussed by numerical simulations in detail.

## Acknowledgments

This work was partially supported by National Natural Science Foundation of
China under grants 62273240, the Natural Science Foundation of Shaanxi
Province (Grant No. 2023-JC-QN-0003), the Fundamental Research Funds for the
Central Universities(Grant No. 23GH020210), Suzhou Municipal Science and
Technology Bureau under Grant (ZXL2023177), and the basic Research Program of
Taicang, Grant (TC2023JC13).

## Conflict of interest

The authors declare no potential conflict of interests.

## Declaration of generative AI and AI-assisted technologies in the writing
process

During the preparation of this work the authors used ChatGPT in order to
improve language and readability. After using this tool, the authors reviewed
and edited the content as needed and take full responsibility for the content
of the publication.

##  Appendix A Proof of proposition1

Proof: The proposition1 can be proven by Lyapunov method by two steps. This
first step is prestabilizing the inject estimation error z. The second step is
ensuring the dynamic scaling $r$ is bounded.

The Lyapunov function candidate about the z is chosen as

$V_{z}=\frac{1}{2}\textbf{{z}}^{T}\textbf{{J}}\textbf{{z}}.$ (36)

The dynamics of $\boldsymbol{\omega}_{EI}^{B}$ can be obtained from (12),
(13a), (14) that

$\dot{\boldsymbol{\omega}}_{EI}^{B}=\textbf{{J}}^{-1}(-\boldsymbol{\omega}_{EI}^{B}\times\textbf{{J}}\boldsymbol{\omega}_{EI}^{B}+\boldsymbol{\tau}_{c}^{B})+4\textbf{{J}}^{-1}\beta(\underline{\varpi})\dot{\boldsymbol{\sigma}}_{2BE}.$
(37)

One can derive the dynamics of $\boldsymbol{\omega}_{BE}^{B}$ from (1b), (10b)
, and (37) that

$\dot{\boldsymbol{\omega}}_{BE}^{B}=\textbf{{J}}^{-1}(\boldsymbol{\omega}_{EI}^{B}\times\textbf{{J}}\boldsymbol{\omega}_{EI}^{B}-\boldsymbol{\omega}_{BI}^{B}\times\textbf{{J}}\boldsymbol{\omega}_{BI}^{B})-4\textbf{{J}}^{-1}\beta(\underline{\varpi})\dot{\boldsymbol{\sigma}}_{2BE}.$
(38)

Then, by invoking (5), (16) (17a), and (17b) and using the inequalities
$\|\boldsymbol{\omega}_{EI}^{B}\|\leq\varpi$,
$\|\textbf{{a}}\|\leq\|\textbf{{a}}-\textbf{{b}}\|+\|\textbf{{b}}\|$, and
$2\textbf{{a}}\textbf{{b}}\leq\|\textbf{{a}}\|^{2}+\|\textbf{{b}}\|^{2}$, the
time derivative of $V_{z}$ along (15) and (38) can be obtained as

$\begin{split}\dot{V}_{z}=&\textbf{{z}}^{T}\textbf{{J}}r^{-1}\textbf{{J}}^{-1}(\boldsymbol{\omega}_{EI}^{B}\times\textbf{{J}}\boldsymbol{\omega}_{EI}^{B}-\boldsymbol{\omega}_{BI}^{B}\times\textbf{{J}}\boldsymbol{\omega}_{BI}^{B})\\\
&-\textbf{{z}}^{T}\textbf{{J}}\left\\{r^{-1}4\textbf{{J}}^{-1}\beta(\underline{\varpi})\dot{\boldsymbol{\sigma}}_{2BE}-\dot{r}^{-1}\boldsymbol{\omega}_{BE}^{B}\right\\}\\\
\leq&\textbf{{z}}^{T}\left(\boldsymbol{\omega}_{EI}^{B}\times\textbf{{J}}\textbf{{z}}-\underline{\beta}(\underline{\varpi})\textbf{{z}}\right)-r^{-1}\textbf{{z}}^{T}\textbf{{J}}\dot{r}\textbf{{z}}\\\
\leq&J_{M}\|\underline{\varpi}\|\|\textbf{{z}}\|^{2}-\underline{\beta}(\underline{\varpi})\|\textbf{{z}}\|^{2}+\dfrac{J_{m}k_{r}}{J_{M}}\|\textbf{{z}}\|^{2}\\\
=&-(1+\rho_{\varpi})\|\textbf{{z}}\|^{2}\end{split}$ (39)

which implies z converges to zeros exponentially.

To show the boundedness of r, consider

$V_{o}=V_{z}+V_{\varpi}+V_{\sigma}+V_{r}$ (40)

where

$V_{\varpi}=\frac{1}{2}(\varpi-\underline{\varpi})^{2}$ $V_{\sigma}=2{\rm
ln}(1+\sigma^{2}_{BE})$ $V_{r}=\dfrac{J_{m}}{2}(r-1)^{2}.$

Using (11) and (37), the dynamics of $\varpi$ is derived as

$\displaystyle\dot{\varpi}=$
$\displaystyle\varpi^{-1}(\boldsymbol{\omega}_{EI}^{B})^{T}\dot{\boldsymbol{\omega}}_{EI}^{B}$
(41) $\displaystyle=$
$\displaystyle\varpi^{-1}(\boldsymbol{\omega}_{EI}^{B})^{T}\textbf{{J}}^{-1}\\!\left\\{\\!-\boldsymbol{\omega}_{EI}^{B}\\!\times\\!\textbf{{J}}\boldsymbol{\omega}_{EI}^{B}+\boldsymbol{\tau}_{c}^{B}+4\beta(\underline{\varpi})\dot{\boldsymbol{\sigma}}_{2BE}\right\\}\\!.$

From (13b) and (41), it follows that

$\dot{\varpi}-\dot{\underline{\varpi}}=\varpi^{-1}(\boldsymbol{\omega}_{EI}^{B})^{T}4\textbf{{J}}^{-1}\beta(\underline{\varpi})\dot{\boldsymbol{\sigma}}_{2BE}-K_{\underline{\varpi}}(\varpi-\underline{\varpi})$
(42)

Take the time derivative of $V_{\varpi}$ along (42), one can obtain

$\displaystyle\dot{V}_{\varpi}=$
$\displaystyle(\varpi\\!-\\!\underline{\varpi})\\!\left\\{\varpi^{-1}(\boldsymbol{\omega}_{EI}^{B})^{T}\\!4\textbf{{J}}^{-1}\beta(\underline{\varpi})\dot{\boldsymbol{\sigma}}_{2BE}\\!-\\!K_{\underline{\varpi}}(\varpi\\!-\\!\underline{\varpi})\right\\}$
(43) $\displaystyle\leq$
$\displaystyle\frac{1}{2}\textbf{{z}}^{2}-(\varpi-\underline{\varpi})^{2}\left\\{K_{\underline{\varpi}}-8\left(\dfrac{\|\boldsymbol{\omega}_{EI}^{B}\|\beta(\underline{\varpi})r}{J_{m}\varpi}\right)^{2}\right\\}$

Applying (14) and (4), one has

$\displaystyle\dot{V}_{\sigma}=$
$\displaystyle\dfrac{4\boldsymbol{\sigma}^{T}_{BE}}{1+\boldsymbol{\sigma}^{2}_{BE}}G(\boldsymbol{\sigma}_{BE})(\boldsymbol{\omega}_{BE}^{B}-K_{\sigma}\boldsymbol{\sigma}_{BE})$
(44) $\displaystyle\leq$
$\displaystyle\frac{1}{2}\textbf{{z}}^{2}-\left(K_{\sigma}-\frac{1}{2}r^{2}\right)\|\boldsymbol{\sigma}_{BE}\|^{2}$

Additionally, following the calculations in (16), one can obtain

$\displaystyle\dot{V}_{r}=$ $\displaystyle
J_{m}(r-1)\left\\{\dfrac{\textit{r}}{J_{m}}(J_{M}\|\varpi-\underline{\varpi}\|)-\dfrac{k_{r}}{J_{M}}(r-1)\right\\}$
(45) $\displaystyle\leq$
$\displaystyle\frac{1}{2}\textbf{{r}}^{2}J_{M}\|\varpi-\underline{\varpi}\|^{2}-\left(\dfrac{J_{m}k_{r}}{J_{M}}-\frac{1}{2}J_{M}\right)(r-1)^{2}$

Finally, differentiating $V_{o}$ along (39), (43), (44), and (45), yields

$\dot{V}_{o}\leq-\rho_{\varpi}\|\textbf{{z}}\|^{2}-\rho_{\underline{\varpi}}(\varpi-\underline{\varpi})^{2}-\rho_{\sigma}\|\boldsymbol{\sigma}_{BE}\|^{2}-\rho_{r}(r-1)^{2}$
(46)

which implies that the system is exponentially stable and $r$ is bounded. This
completes the proof.

## References

  * [1] C. Wei, Q. Chen, J. Liu, Z. Yin, J. Luo, An overview of prescribed performance control and its application to spacecraft attitude system, Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering 235 (4) (2021) 435–447.
  * [2] R. Qi, X. Dong, D. Chao, Y. Wang, Constrained attitude tracking control and active sloshing suppression for liquid-filled spacecraft, ISA Transactions 132 (2023) 292–308.
  * [3] S. Nakano, T. W. Nguyen, E. Garone, T. Ibuki, M. Sampei, Explicit reference governor on so(3) for torque and pointing constraint management, Automatica 155 (2023) 111103.
  * [4] Y. Qu, X. Zhong, F. Zhang, X. Tong, L. Fan, L. Dai, Robust disturbance observer-based fast maneuver method for attitude control of optical remote sensing satellites, Acta Astronautica 201 (2022) 83–93.
  * [5] H. B. Hablani, Attitude commands avoiding bright objects and maintaining communication with ground station, Journal of Guidance, Control, and Dynamics 22 (6) (1999) 759–767.
  * [6] H. Cui, X. Cheng, Anti-unwinding attitude maneuver control of spacecraft considering bounded disturbance and input saturation, SCIENCE CHINA Technological Sciences 55 (9) (2012) 2518–2529.
  * [7] Q. Dang, K. Liu, J. Wei, Explicit reference governor based spacecraft attitude reorientation control with constraints and disturbances, Acta Astronautica 190 (2022) 455–464.
  * [8] S. H. Pourtakdoust, M. F. Mehrjardi, M. Hajkarim, Attitude estimation and control based on modified unscented kalman filter for gyro-less satellite with faulty sensors, Acta Astronautica 191 (2022) 134–147.
  * [9] J. Hu, H. Zhang, Bounded output feedback of rigid-body attitude via angular velocity observers, Journal of Guidance, Control, and Dynamics 36 (4) (2013) 1240–1248.
  * [10] H. C. Kjellberg, E. G. Lightsey, Discretized constrained attitude pathfinding and control for satellites, Journal of Guidance, Control, and Dynamics 36 (5) (2013) 1301–1309.
  * [11] B. Wie, D. Bailey, C. Heiberg, Rapid multitarget acquisition and pointing control of agile spacecraft, Journal of Guidance Control Dynamics 25 (1) (2002) 96–104.
  * [12] Y. Cheng, D. Ye, Z. Sun, S. Zhang, Spacecraft reorientation control in presence of attitude constraint considering input saturation and stochastic disturbance, Acta Astronautica 144 (2018) 61–68.
  * [13] S. Kulumani, T. Lee, Constrained geometric attitude control on so (3), International Journal of Control, Automation and Systems 15 (6) (2017) 2796–2809.
  * [14] B. Hua, J. He, H. Zhang, Y. Wu, Z. Chen, Spacecraft attitude reorientation control method based on potential function under complex constraints, Aerospace Science and Technology 144 (2024) 108738.
  * [15] U. Lee, M. Mesbahi, Feedback control for spacecraft reorientation under attitude constraints via convex potentials, IEEE Transactions on Aerospace and Electronic Systems 50 (4) (2014) 2578–2592.
  * [16] Q. Hu, B. Chi, M. R. Akella, Anti-unwinding attitude control of spacecraft with forbidden pointing constraints, Journal of Guidance, Control, and Dynamics 42 (4) (2019) 822–835.
  * [17] Q. Shen, C. Yue, C. H. Goh, B. Wu, D. Wang, Rigid-body attitude stabilization with attitude and angular rate constraints, Automatica 90 (2018) 157 – 163.
  * [18] U. V. Kalabić, R. Gupta, S. D. Cairano, A. M. Bloch, I. V. Kolmanovsky, MPC on manifolds with an application to the control of spacecraft attitude on SO(3), Automatica 76 (2017) 293 – 300.
  * [19] D. Y. Lee, R. Gupta, U. V. Kalabić, S. Di Cairano, A. M. Bloch, J. W. Cutler, I. V. Kolmanovsky, Geometric mechanics based nonlinear model predictive spacecraft attitude control with reaction wheels, Journal of Guidance, Control, and Dynamics 40 (2) (2017) 309–319.
  * [20] M. M. Nicotra, E. Garone, The explicit reference governor: A general framework for the closed-form control of constrained nonlinear systems, IEEE Control Systems Magazine 38 (4) (2018) 89–107.
  * [21] E. Garone, M. M. Nicotra, Explicit reference governor for constrained nonlinear systems, IEEE Transactions on Automatic Control 61 (5) (2016) 1379–1384.
  * [22] M. M. Nicotra, D. Liao-McPherson, L. Burlion, I. V. Kolmanovsky, Spacecraft attitude control with nonconvex constraints: An explicit reference governor approach, IEEE Transactions on Automatic Control 65 (8) (2020) 3677–3684.
  * [23] M. Hosseinzadeh, E. Garone, An explicit reference governor for the intersection of concave constraints, IEEE Transactions on Automatic Control 65 (1) (2020) 1–11.
  * [24] M. Nicotra, E. Garone, Constrained control of nonlinear systems: The explicit reference governor and its application to unmanned aerial vehicles, Ph.D. thesis, Universite Libre de Bruxelles (2016).
  * [25] B. Chi, Q. Hu, Saturated explicit reference governor for spacecraft constrained attitude reorientation control, Aerospace Science and Technology 145 (2024) 108874\.
  * [26] M. N. Hasan, Y. Chen, J. Liang, A. Wen, Fixed-time fault-tolerant attitude control for flexible spacecraft without angular velocity sensor, ISA Transactions (2023, in press).
  * [27] S. Yang, M. R. Akella, F. Mazenc, Immersion and invariance observers for gyro-free attitude control systems, Journal of Guidance Control Dynamics 39 (11) (2016) 2567–2574.
  * [28] E. Espíndola, Y. Tang, A new angular velocity observer for attitude tracking of spacecraft, ISA Transactions 130 (2022) 377–388.
  * [29] R. Schlanbusch, E. I. Ingar Grotli, Hybrid certainty equivalence control of rigid bodies with quaternion measurements, IEEE Transactions on Automatic Control 60 (9) (2015) 2512–2517.
  * [30] J. G. Romero, R. Ortega, I. Sarras, A globally exponentially stable tracking controller for mechanical systems using position feedback, IEEE Transactions on Automatic Control 60 (3) (2015) 818–823.
  * [31] Q. Dang, H. Gui, H. Wen, Dual-quaternion-based spacecraft pose tracking with a global exponential velocity observer, Journal of Guidance, Control, and Dynamics 42 (9) (2019) 2106–2115.
  * [32] Q. Shen, C. Yue, C. H. Goh, Velocity-free attitude reorientation of a flexible spacecraft with attitude constraints, Journal of Guidance, Control, and Dynamics 40 (5) (2017) 1293–1299.
  * [33] J. L. Junkins, H. Schaub, Analytical mechanics of space systems, American Institute of Aeronautics and Astronautics, 2009.
  * [34] H. Gui, G. Vukovich, Adaptive integral sliding mode control for spacecraft attitude tracking with actuator uncertainty, Journal of the Franklin Institute 352 (12) (2015) 5832 – 5852.
  * [35] M. Diaz Ramos, H. Schaub, Kinematic steering law for conically constrained torque-limited spacecraft attitude control, Journal of Guidance, Control, and Dynamics 41 (9) (2018) 1990–2001.
  * [36] A. Bemporad, Reference governor for constrained nonlinear systems, IEEE Transactions on Automatic Control 43 (3) (1998) 415–419.


11institutetext: Department of Computer Science, Aalborg University, Aalborg,
Denmark 11email<EMAIL_ADDRESS>
22institutetext: Unit of Business Intelligence, North Denmark Region, Aalborg,
Denmark 22email<EMAIL_ADDRESS>33institutetext: TU Wien, Vienna, Austria
33email<EMAIL_ADDRESS>
# Hospitalization Length of Stay Prediction using Patient Event Sequences

Emil Riis Hansen(✉) 11 0000-0003-4103-1244 Thomas Dyhre Nielsen 11
0000-0002-4823-6341 Thomas Mulvad 22 0000-0002-4457-3965 Mads Nibe Strausholm
22 0000-0002-4099-3209 Tomer Sagi 11 0000-0002-8916-0128 Katja Hose 1133
0000-0001-7025-8099

###### Abstract

Predicting patients’ hospital length of stay (LOS) is essential for improving
resource allocation and supporting decision-making in healthcare
organizations. This paper proposes a novel approach for predicting LOS by
modeling patient information as sequences of events. Specifically, we present
a transformer-based model, termed Medic-BERT (M-BERT), for LOS prediction
using the unique features describing patients’ medical event sequences. We
performed empirical experiments on a cohort of more than $45k$ emergency care
patients from a large Danish hospital. Experimental results show that M-BERT
can achieve high accuracy on a variety of LOS problems and outperforms
traditional non-sequence-based machine learning approaches.

###### Keywords:

length of stay prediction transformers sequence models

##  1 Introduction

Increasingly scarce hospital resources challenge (often oversaturated)
hospital wards, with a negative impact on the quality of health care at the
hospitals [1] Models for predicting the remaining time of patient admissions,
i.e., patient length of stay (LOS), could be invaluable for healthcare
facilities to plan the availability of beds, staff, and other essential
resources. For instance, automatic prediction of discharge time could be used
in administrative planning systems for preemptively freeing in-hospital
resources to alleviate hospital ward oversaturation [20]. However, LOS
prediction is a challenging problem, requiring methods for handling missing
data [18] and integration of temporal event dependencies.

Previous work on LOS prediction models patient hospitalizations using tabular
data with imputation techniques for replacing missing values [2]. While
tabular data is the most common data representation in Machine Learning (ML)
models, it has several drawbacks. Among others, it does not provide immediate
support for integrating the temporal dependencies between observations, such
as the order of treatments, or the time of conducted procedures. Moreover,
standard ML techniques for tabular data, such as artificial neural networks
(ANN), gradient boosting (GB), and support vector machines (SVMs), require
complete data, hence often relying on imputation techniques when data is
incomplete. However, missing data observations in healthcare data are often
not missing at random (NMAR), meaning that the mere fact that an observation
is missing is in itself important information [11].

To alleviate the problem of temporal dependencies and missing data, attention
models, also known as transformers, have recently been investigated for
Electronic Health Record (EHR) data formatted as sequences of medical events
[10]. However, embedding-based transformer approaches have, to the best of our
knowledge, not previously been applied to LOS prediction. This work examines
how attention models, can be utilized for this task.

Attention models using self-attention alleviate the inefficiency of recurrence
networks for long sequences [19]. However, they still capture significant
sequential information by learning from the order of tokens in the sequence.
In medical data, multiple observations may be given the same timestamp, with
no meaning assigned to their individual order within the corresponding event.
For example, a blood panel drawn from a patient contains several individual
measurements whose internal order is insignificant. Based on layers of
transformer encoders, we propose a revised attention model, henceforth called
Medic-BERT (M-BERT), based on the Natural Language Processing (NLP) model BERT
[6] and its revised semi-supervised training method. We employ the model for
LOS prediction based on sequences of patient specific medical events happening
during hospitalization and which exhibit the event concurrences common in
patient data. We evaluate our method on a cohort of more than $45k$ patient
admissions from a large Danish hospital with diverse medical events, such as
vital measurements, medication administration, laboratory tests, and conducted
procedures.

The rest of this paper is structured as follows. Section 2 reviews related
work on LOS prediction. Section 3 presents our model for representing
hospitalizations as event sequences and our proposed M-BERT model for this
unique data. Section 4 describes the large dataset used to evaluate this work
and Section 5 the evaluation and its results. We conclude with Section 6. This
work is an extended version of Hansen et al. [7] published at AIME 2023.

##  2 Related Work

The LOS prediction problem can be stated in different ways varying in the
required resolution of the prediction, applied methods, and structure of
patient data. Iwase et al. [9] investigate the use of Machine Learning (ML)
methods for LOS prediction of intensive care unit patients (ICU). Using Random
Forest (RF), Gradient Boosting (GB), and Artificial Neural Networks (ANN)
technologies, on binary stratification of patient admissions into long (more
than one week) and short (less than one week) stays. RF showed high predictive
performance. In contrast to the tabular structure of their data, we instead
attempt to utilize the temporal dependencies within the data by structuring
the patient data as medical event sequences. In Batista et al. [4], patients
are stratified into three categories (LOS $<3$, $3<$ LOS $<10$, LOS $>=10$)
using RF and Support Vector Machines (SVM). The highest performance was
achieved for RF models. While the nature of their ML models requires manual
feature selection and missing data imputation, our method does not require any
feature selection for in-hospital measurements, nor does it rely on data
imputation.

The most challenging LOS prediction task is regression-based LOS, where the
precise hospitalization LOS as a real value is to be predicted. Barsasella et
al. [3] investigated a range of classical AI models, such as Decision Trees
(DTs), RFs, Logistic Regression (LR), and SVMs for real-valued LOS prediction.
While standard ML models work for tabular structured patient data, they
neglect the important temporal dependencies between in-hospital patient
events. While much work has been done in LOS prediction and regression using
standard ML methods, only a few works attempt the usage of sequence models.

Seminal work in the area of transformer models for patient sequence data
includes that of Li et al. [12] and Rasmy et al. [16], where sequences of
diagnosis codes for consecutive hospitalizations are used as input to a BERT
like transformer model for diagnosis code prediction tasks. While Rasmy et al.
[16] investigate prolonged LOS $>7$ days event prediction as a model pre-
training task, they do not predict future LOS for new patients. Furthermore,
while their approaches are only investigated for sequences consisting of
diagnosis codes, we evaluate sequences consisting of multiple patient event
types.

In Meng el al. [13], a transformer model is used to predict the onset of
depression. They integrate various medical types for patient event sequences,
including diagnosis and medication events. While their sequences include code-
based medical events, we extend with medical measurement events, such as
laboratory tests and vital measures while encoding the event value as part of
the event token.

The work closest to ours is that of Song et al. [19], where an attention model
is used to predict categorical LOS based on patient sequences. Whereas their
approach targets dense time-series data, where the same measurements are
present at every timestep, instead, we investigate a learned embedding
approach where the input sequences consist of codes with different types of
medical events.

##  3 Transformer Models for EHR

In this section, we describe how EHR data can be modeled as medical event
sequences and our usage of transformer models for learning from these
sequences.

###  3.1 Hospitalizations as Event Sequences

Patient hospitalizations can be naturally modeled as sequences of medical
procedures for determining, measuring, or diagnosing the patient’s condition.
Other medical procedures are therapeutic and intend to treat or cure the
patient. To standardize how medical procedures are described, medical
facilities code procedure concepts using accepted medical taxonomies often
used in ML applications [8], such as the Anatomical Therapeutic Classification
(ATC) [17] for medication administration and the International Classification
of Diseases and Related Health Problems (ICD) [5] for condition diagnosis.
Hence, a patient hospitalization can be described as the sequence of concept
tokens detailing medical procedures pertaining to a patient coded using
medical concepts from accepted medical taxonomies. An example patient
hospitalization sequence is illustrated in Figure 1.

Furthermore, as a patient’s medical history is crucial for correct management
and treatment, we pre-pend the patient’s medical history to the
hospitalization sequence as a tokenized vector as illustrated in Figure 1. The
vector consists of 38 tokens describing the patient’s medical history,
including comorbidities from the Charlson Index [21], five years of
medications prescription history grouped by the first level of the ATC
hierarchy [17], and the mode, time, and initial triage category [23] of
hospitalization. The historical information included in the vector is
summarized in Table 1.

Figure 1: Illustration of a patient event sequence. Starting with the patient’s medical history as summarized in Table 1, the patient is initially diagnosed with the ICD-10 code _Z039_. Subsequent vital measurements and laboratory tests are performed to monitor the patient’s state and determine the patient’s underlying condition. Consequently, the patient is diagnosed with acute cystitis without hematuria (ICD-10 code _N300_), and antibiotic treatment is initiated with nitrofurantoin (ATC code _J01XE01_). After additional procedure and therapeutic medical events, the patient is released from the hospital. Table 1: Patient medical history. Data |  #Tokens
---|---
Comorbidities |  18
Prescription history |  14
Mode, time, & triage |  6

###  3.2 Hospital Measurement Events

For some medical procedures, such as vital measurements and laboratory tests,
a numerical measurement value accompanies the procedure. While other works in
transformer models for EHR data disregard the numerical values of measurements
[13, 16, 12], we instead integrate this information as part of the patient
input sequences because numerical measurement values add important information
regarding the state of a patient. For example, the knowledge that a
temperature measurement was performed is naturally important information.
Still, from the measurement value of $40.1^{\circ}C$, we learn that the
patient has a fever. Using patient-specific threshold values for measurements
based on age, gender, and pregnancy status, we map measurement values into
tokens representing either normal, abnormal-low, or abnormal-high findings.
For example, given that we measure an albumin level of $56$ g/L for a
$31$-year-old male patient, we would create the token _albumin-high_ to
reflect that the value of the measurement was above what is expected
($36$-$48$ g/L) for a patient with the given demography.

Figure 2: Medical event sequence pre-pended with the patient’s medical
history.

###  3.3 Transformer Models for EHR Data

To investigate the challenge of LOS, we examine a modified version, henceforth
termed Medic-BERT (M-BERT), of the Bidirectional Encoder Representations from
Transformers (BERT) [6] model for EHR data. BERT is an NLP model based on a
stack of encoder layers from the transformer architecture introduced by
Vaswani et al. [22]. We argue that sequence models, such as BERT exhibit
properties beneficial for solving medical tasks based on EHR data. The
transformer encoder naturally handles the complex long-term dependencies that
occur between medical concepts through its utilization of multi-head self-
attention. BERT can naturally integrate disparate modalities, such as
diagnostic end therapeutic events, as each event is encoded as an
n-dimensional vector token. Furthermore, BERT naturally operates in domains
with irregular intervals between events, as is the case with EHR data. We,
therefore, investigate our modified version of BERT for LOS prediction for
patient event sequences. The M-BERT architecture is illustrated in Figure 2.

M-BERT learns an embedding for each medical event token while trained toward
LOS prediction. The position embedding enables the model to learn from the
temporal dependencies within a sequence. We use a static position embedding as
described in Waswani et al. [22]; however, modified for the usage on medical
event sequences as described in Section 3.4. As patient demographics are a
vital part of any medical prediction model, we pay special attention to this
information by adding a trainable age and sex embedding at each medical event
[12]. Furthermore, as in the original BERT model, we use the special
classification (CLS) token as a final aggregate representation of the
sequence. Hence, as illustrated in Figure 2, the CLS representation is fed to
a linear output layer for LOS classification and regression tasks.

###  3.4 Position Embedding

Due to the nature of patient care and hospital administration, some
measurement events tend to chunk together. Clinicians will often conduct
various measurements, such as blood pressure, temperature, and heart rate,
over a short time and later persist the information into the patient’s EHR.
Hence, we are sometimes prevented from knowing the specific times and order of
medical events. This effect is most frequent for vital measurements and
laboratory tests. Multiple laboratory tests are often conducted on the same
patient sample, such as a single blood sample, making it impossible to know
the chronological order for such measurement events. An example sequence of
patient events chunking together is illustrated in Figure 3. To enable the
model to understand that these events have no fixed ordering, we assign the
same position embedding to events co-located in time as illustrated in the
position embedding of Figure 2. For example, the $L_{1}$ and $L_{2}$ events
are both mapped to the $E_{P_{3}}$ position embedding.

Figure 3: Patient event sequence illustrating events grouping together.

##  4 Data

We compiled a Danish dataset of patients admitted to a large hospital in
northern Jutland from the period 2018-2021 to investigate transformer models
for LOS prediction. The dataset consists of $48,177$ emergency care patients
with admissions longer than one day. Figure 4 illustrates the distribution of
the remaining length of stay for the patient cohort. In a clinical setting of
resource allocation for emergency care patients, we are mostly interested in
the planning of patient care for a couple of weeks in advance. Hence, we clip
long patient stays to 30 days of admission to better fit and optimize for the
clinical setting, as we are not interested in precisely predicting the long-
tail distribution of the data.

In Denmark, each person can be uniquely identified by an identification number
from the central person register (CPR) henceforth referred to as a CPR-number.
As medical events are each associated with a unique CPR-number, it is possible
to connect the medical data pertaining to a patient from disparate databases.
As mentioned in Section 3.1, we divide patient medical data into historic
information and admission-specific information. Historic information includes
prescriptions, comorbidities, and mode and time of hospitalization, whereas
admission information includes laboratory tests, vital measurements, hospital-
administered medicine, and procedure codes.

(a) Histogram of length of stay.

(b) Categorical length of stay.

Figure 4: Illustrations for length of stay distribution over patient hospitalizations. Table 2: Concept types with their occurrences in the dataset. Event Type |  Tokens | Occurrences
---|---|---
Lab Tests |  $748$ | $2,774,790$
Vital Measures |  $22$ | $837,931$
Medication |  $1,441$ | $376,591$
Procedures |  $2049$ | $247,924$
History |  $81$ | $1,880,580$

Danish laboratory tests are coded using the Nomenclature for Properties and
Units (NPU) terminology [15]. NPU ensures that laboratory tests are
standardized and patient examinations can be used and understood by all
clinicians. As detailed in Section 3.2, we integrate the semantics of
laboratory test results as part of the event token by appending either high,
low, or normal with respect to the patient demographics to the laboratory
event token, thus encoding the meaning of the result. Vital measurements are
extracted from a system called Clinical Suite and consist of the 7 most common
vital observations, including temperature, oxygen saturation, BMI, pulse,
respiration rate, blood pressure (systolic and diastolic), and oxygen
supplement. As with laboratory tests, patient-specific thresholds are used to
encode the meaning of the result into event tokens. In-hospital administered
medication is coded using the ATC taxonomy [17] and consists of more than $5$k
chemical substances. Procedure codes specify in-hospital procedures performed
on patients. While diagnosis events are an important medical modality, the
dataset does not contain the time of such events. Hence, we omit diagnosis
events from patient sequences. The event types with distinct tokens and total
occurrence in the data are summarized in Table 2.

##  5 Empirical Evaluation and Results

In this section, we explain the experimental settings and results of the
empirical evaluation.

###  5.1 Experimental Setting

Transformer models are trained, as a rule, using unsupervised pre-training for
learning general token embeddings, followed by supervised fine-tuning
targeting a specific downstream task. As we are only focusing on a single
prediction task, we directly train model parameters and token embeddings
toward the downstream task of LOS prediction without pre-training. The
experimental code is available online111https://github.com/dkw-
aau/medic_transformer.

We use the admission data gathered within the first 24 hours of admissions for
LOS predictions. Hence, we remove sequence events happening after 24 hours of
admission. We evaluate our approach on three LOS experiments of increasing
complexity. The first two experiments are a Binary classification of LOS $>2$
days, and a three-class Category task of LOS $>2$, $2\leq$ LOS $\leq 7$, and
LOS $>7$ days with class balances as illustrated in Figure 4b. The last
experiment, termed Real, is a regression task with the objective of predicting
the LOS as the real number, with a histogram of admission times as illustrated
in Figure 4a.

We evaluate our approach against three ML models: RF, ANN, and SVM. We use
implementations of the models from the Sklearn library [14] with standard
model hyperparameters. In preparing samples for these models, we use the
latest measured value for each event type (within 24 hours of admission) as
input features [9]. Subsequently, the mean of variables is used for imputation
of missing data values, and variables are scaled to values between 0 and 1.
Lastly, a chi2 test is used in feature selection for the selection of the 50
most relevant features.

Our model is trained on a random split distributed as 80/10/10 of all patient
samples for training, validation, and testing. We use the loss of the
evaluation data for early stopping training if the loss does not go down
within ten training epochs. The model architecture has six hidden layers with
an intermediate layer size of 288, eight attention heads, and input token
embeddings have a size of 288. We truncate sequences to 256 tokens, as most
sequences adhere to this limit. To counter overfitting, we add a dropout layer
with a probability of 10% after the output of the final encoder layer,
attention dropout at every layer with a probability of 10%, and weight decay
of $0,003$. Furthermore, experiments were performed with a learning rate of
1e-5.

(a) AUROC curves for the Binary LOS experiment over all methods.

(b) Binary experiment AUROC curves for different age groups over M-BERT.

Figure 5: AUROC plots for experimental results on the binary prediction task.

###  5.2 Results

Table 3 presents the Area Under the Receiver Operating Characteristics (AUROC)
and harmonic mean of precision and recall (F1) values for the Binary and
Category experimental settings and Mean Absolute Error (MAE) and Mean Squared
Error (MSE) for the LOS regression task. Furthermore, Figure 5a presents the
AUROC curves for the Binary experimental setting. M-BERT outperforms the
traditional ML techniques in all experimental settings. The results indicate
that transformer models may be able to leverage the temporal dependencies
inherent to patient EHR data for increased predictive accuracy.

Table 3: Experimental results. | Binary |  | Category |  | Real |
---|---|---|---|---|---|---
| AUROC | F1 |  | AUROC | F1 |  | MAE | MSE |
RF | $0.72$ | $0.70$ |  | $0.66$ | $0.45$ |  | $4.18$ | $39.08$ |
ANN | $0.67$ | $0.68$ |  | $0.63$ | $0.43$ |  | $4.09$ | $38.10$ |
SVM | $0.70$ | $0.70$ |  | $0.65$ | $0.38$ |  | $3.56$ | $43.36$ |
M-BERT | 0.78 | 0.77 |  | 0.74 | 0.54 |  | 3.42 | 37.48 |

Furthermore, being a transformer-based model, M-BERT overcomes the challenge
of missing data and imputation, as patient sequences are not required to
contain the same medical events, nor the same sequence length. Figure 5b
presents the AUROC for binary LOS prediction stratified by age groups for
those with more than 100 patient samples. Interestingly, the model performance
is stable for different age groups, with an AUROC of 0.78 for the age groups
60-70 and 70-80 years. Furthermore, stratification based on the sex of
patients yielded similar results, pointing to the robustness of M-BERT for LOS
classification.

##  6 Conclusion

We have presented a novel approach for predicting LOS by modeling patient
information as event sequences. Our approach adapts the transformer machine-
learning approach for sequence prediction, which is able to handle the unique
features of medical event sequences, namely grouped events and a variety of
data types. Our empirical evaluation on a large cohort of emergency care
patients from a Danish hospital demonstrates that our model can achieve high
accuracy on various LOS problems, while outperforming traditional non-sequence
machine learning approaches. Future work could include pre-training of the
transformer-based model on a medical task to further improve its performance.
Overall, the proposed approach has the potential to improve resource
allocation and support decision making in healthcare organizations by
providing accurate predictions of LOS. All experimental code is available
online.

## Acknowledgments

This work is partially supported by the Poul Due Jensen Foundation and the
Region North Denmark Health Innovation Foundation.

## References

  * [1] af Ugglas, B., Djärv, T., Ljungman, P.L., Holzmann, M.J.: Association between hospital bed occupancy and outcomes in emergency care: a cohort study in stockholm region, sweden, 2012 to 2016. Ann. Emerg. Med. 76(2), 179–190 (2020)
  * [2] Bacchi, S., Tan, Y., Oakden-Rayner, L., Jannes, J., Kleinig, T., Koblar, S.: Machine learning in the prediction of medical inpatient length of stay. Intern. Med. J. 52(2), 176–185 (2022)
  * [3] Barsasella, D., Bah, K., Mishra, P., Uddin, M., Dhar, E., Suryani, D.L., et al.: A machine learning model to predict length of stay and mortality among diabetes and hypertension inpatients. Medicina 58(11), 1568 (2022)
  * [4] Batista, R.W., Sanchez-Arias, R.: A methodology for estimating hospital intensive care unit length of stay using novel machine learning tools. In: ICMLA ’20. pp. 827–832. IEEE (2020)
  * [5] Cartwright, D.J.: Icd-9-cm to icd-10-cm codes: what? why? how? (2013)
  * [6] Devlin, J., Chang, M., Lee, K., Toutanova, K.: BERT: pre-training of deep bidirectional transformers for language understanding. In: NAACL-HLT ’19. pp. 4171–4186 (2019). https://doi.org/10.18653/v1/n19-1423
  * [7] Hansen, E.R., Nielsen, T.D., Mulvad, T., Strausholm, M.N., Sagi, T., Hose, K.: Patient event sequences for predicting hospitalization length of stay. In: Proceedings of the 21st International Conference on Artificial Intelligence in Medicine (2023)
  * [8] Hansen, E.R., Sagi, T., Hose, K., Lip, G.Y.H., Larsen, T.B., Skjøth, F.: Assigning diagnosis codes using medication history. Artif. Intell. Medicine 128, 102307 (2022)
  * [9] Iwase, S., Nakada, T.a., Shimada, T., Oami, T., Shimazui, T., Takahashi, N., et al.: Prediction algorithm for icu mortality and length of stay using machine learning. Scientific reports 12(1), 1–9 (2022)
  * [10] Lequertier, V., Wang, T., Fondrevelle, J., Augusto, V., Duclos, A.: Hospital length of stay prediction methods: A systematic review. Med. Care 59(10), 929–938 (2021)
  * [11] Li, J., Yan, X.S., Chaudhary, D., Avula, V., Mudiganti, S., Husby, H., et al.: Imputation of missing values for electronic health record laboratory data. NPJ digital medicine 4(1), 1–14 (2021)
  * [12] Li, Y., Rao, S., Solares, J.R.A., Hassaine, A., Ramakrishnan, R., Canoy, D., Zhu, Y., Rahimi, K., Salimi-Khorshidi, G.: Behrt: transformer for electronic health records. Scientific reports 10(1), 1–12 (2020)
  * [13] Meng, Y., Speier, W., Ong, M.K., Arnold, C.W.: Bidirectional representation learning from transformers using multimodal electronic health record data to predict depression. IEEE J. Biomed. Inform. 25(8), 3121–3129 (2021)
  * [14] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O.e.a.: Scikit-learn: Machine learning in Python. JMLR 12, 2825–2830 (2011)
  * [15] Pontet, F., Petersen, U.M., Fuentes-Arderiu, X., Nordin, G., Bruunshuus, I., Ihalainen, J., et al.: Clinical laboratory sciences data transmission: the npu coding system. Stud. Health Technol. Inform. 150, 265 (2009)
  * [16] Rasmy, L., Xiang, Y., Xie, Z., Tao, C., Zhi, D.: Med-bert: pretrained contextualized embeddings on large-scale structured electronic health records for disease prediction. NPJ digital medicine 4(1), 1–13 (2021)
  * [17] Ronning, M.: A historical overview of the atc/ddd methodology. WHO drug information 16(3), 233 (2002)
  * [18] Sha, Y., Wang, M.D.: Interpretable predictions of clinical outcomes with an attention-based recurrent neural network. In: ACM BCB ’17. pp. 233–240 (2017)
  * [19] Song, H., Rajan, D., Thiagarajan, J., Spanias, A.: Attend and diagnose: Clinical time series analysis using attention models. In: AAAI ’18. vol. 32 (2018)
  * [20] Stone, K., Zwiggelaar, R., Jones, P., Mac Parthaláin, N.: A systematic review of the prediction of hospital length of stay: Towards a unified framework. PLOS Digital Health 1(4), e0000017 (2022)
  * [21] Sundararajan, V., Henderson, T., Perry, C., Muggivan, A., Quan, H., Ghali, W.A.: New icd-10 version of the charlson comorbidity index predicted in-hospital mortality. Journal of clinical epidemiology 57(12), 1288–1294 (2004)
  * [22] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., el al.: Attention is all you need. Adv. Neural Inf. Process. Syst. 30 (2017)
  * [23] Wireklint, S.C., Elmqvist, C., Göransson, K.E.: An updated national survey of triage and triage related work in sweden: a cross-sectional descriptive and comparative study. SJTREM ’21 29(1), 1–8 (2021)