# Ultra Local Nonlinear Unknown Input Observers
for Robust Fault Reconstruction

Farhad Ghanipoor1, Carlos Murguia1, Peyman Mohajerin Esfahani2, and Nathan van
de Wouw1 1 Farhad Ghanipoor, Carlos Murguia, and Nathan van de Wouw are with
the Mechanical Engineering Department, Eindhoven University of Technology,
Eindhoven, The Netherlands. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>and<EMAIL_ADDRESS>Peyman Mohajerin Esfahani is with the Delft Center for
Systems and Control, Delft University of Technology, Delft, The Netherlands.
Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we present a methodology for actuator and sensor fault
estimation in nonlinear systems. The method consists in augmenting the system
dynamics with an approximated ultra-local model (a finite chain of
integrators) for the fault vector and constructing a Nonlinear Unknown Input
Observer (NUIO) for the augmented dynamics. Then, fault reconstruction is
reformulated as a robust state estimation problem in the augmented state (true
state plus fault-related state). We provide sufficient conditions that
guarantee the existence of the observer and stability of the estimation error
dynamics (asymptotic stability of the origin in the absence of faults and ISS
guarantees in the faulty case). Then, we cast the synthesis of observer gains
as a semidefinite program where we minimize the $\mathcal{L}_{2}$-gain from
the model mismatch induced by the approximated fault model to the fault
estimation error. Finally, simulations are given to illustrate the performance
of the proposed methodology.

##  I INTRODUCTION

The increasing expectation of product quality and cost-efficient operation of
engineered systems has led to an increasing demand on process reliability.
This can be achieved through predictive maintenance technology, for which
methods and techniques for fault detection (is a fault occurring?), isolation
(what is the fault source?), and estimation (how large the fault is?) are
fundamental ingredients. Fault Detection and Isolation (FDI) has been an
active field of research for many years (see, for instance, [1]-[8], and
references therein). However, most solutions have been developed for
engineered systems with linear dynamics, when in practice many of them are
highly nonlinear in nature (e.g., robotics, power, transportation, water, and
manufacturing). The standard path to tackle nonlinear behaviour is to obtain
approximated linear models of the system and then apply linear FDI techniques.
Even though this might provide insight about the presence and source of the
fault, using linear models on nonlinear systems often leads to high false
negatives/positives since linear FDI schemes struggle to distinguish between
true faults and model mismatches. The latter has motivated the development of
nonlinear FDI schemes built around different ideas [5]-[8], e.g., observer-
based [5],[6], optimization-based [7], and differential-geometric methods [8].
Most of these results provide elegant FDI solutions; however, they impose
strong hard to verify conditions on the system dynamics – mainly fault-output
decoupling by filtering or observability with respect to different subsets of
outputs – which significantly limits their applicability. Moreover, despite
the great value of existing FDI schemes, what is mainly required for
predictive maintenance is fault estimation. That is, we not only need to know
the presence and source of the fault, but also its severity. If the fault is
small and/or slowly growing in magnitude (slow compared to the system time-
scale), predictions of the fault severity can be performed and predictive
maintenance can be scheduled accordingly. The latter is only possible if we
can reconstruct fault signals (at least their magnitude) with the available
information (inputs, outputs, and system models). Note that by estimating the
fault vector, we automatically have detection and isolation by looking at the
support (the nonzero entries) of the estimated vector.

As it is the case for FDI, fault estimation results are fairly mature for
linear systems (see, e.g., [9, 10] for results on linear stochastic and
switchings systems) and still developing for nonlinear dynamics. In [11], the
authors address the problem for nonlinear systems with uniformly Lipschitz
nonlinearities, process faults only (i.e., no sensor faults), and assume the
so-called matching condition (the rank of the fault distribution matrix is
invariant under left multiplication by the output matrix). For this
configuration, they provide an adaptive filter capable of approximately
reconstructing the actuator fault vector. The matching condition, however, is
a strong assumption, it makes the problem tractable but significantly reduces
the class of system for which results are applicable. In [12], a fault
estimation scheme is introduced for both sensor and process faults using
Nonlinear Unknown Input Observers (NUIO), adaptive Radial Basis Function
Neural Networks (RBFNN), and assuming the matching condition. They prove their
scheme provides boundedness of fault estimation errors.

The authors in [13] do not assume the matching condition. They consider
Lipschitz nonlinearities, simultaneous sensor and process faults, and assume a
standard fault observability condition [14] on the linear part of the
dynamics. They tackle the problem using the notion of intermediate observers,
which consists on having two dedicated observers, one that estimates the fault
and the other the state. Their scheme guarantees bounded fault estimation
errors. In [15], simultaneous additive and multiplicative process faults are
considered. They address the fault estimation problem by decoupling process
nonlinearities and perturbations from the estimation filter dynamics, and
using regression techniques to approximately reconstruct fault signals.
Decoupling nonlinearities leads to linear filters for which linear techniques
can be used to reconstruct fault signals. However, decoupling conditions
impose strong assumptions on the system dynamics, which significantly limits
the applicability of these results.

We remark that the above mentioned results for nonlinear systems guarantee
approximate reconstruction of fault vectors only, i.e., they ensure bounded
estimation errors which, if small enough, still lead to a fairly good estimate
of the true fault. Not having internal models of fault signals makes it
challenging to enforce zero error fault estimation. To this end, we propose a
fault estimation scheme, for process and sensor faults, uniformly Lipschitz
nonlinearities, an without assuming the matching condition, that incorporates
an approximated internal model of the fault vector. We use the notion of
_ultra-local models_ [16]-[17] – which refers to a class of phenomenological
models that are only valid for very short time intervals – to characterize
internal fault models. We then extend the system dynamics with an approximated
internal fault model and construct a NUIO to jointly reconstruct faults and
states. The fault estimation problem is re-formulated as a robust state
estimation problem in the augmented state (true state plus fault-related
state). We provide sufficient conditions that guarantee the existence of the
observer and stability of the estimation error dynamics (asymptotic stability
of the origin in the absence of faults and ISS guarantees in the faulty case).
Then, we cast the synthesis of observer gains as a semidefinite program where
we minimize the $\mathcal{L}_{2}$-gain from the model mismatch induced by the
approximated fault model to the fault estimation error.

Notation: The symbol $\mathbb{R}^{+}$ denotes the set of nonnegative real
numbers. The $n\times n$ identity matrix is denoted by $I_{n}$ or simply $I$
if $n$ is clear from the context. Similarly, $n\times m$ matrices composed of
only zeros are denoted by $\boldsymbol{0}_{n\times m}$ or simply
$\boldsymbol{0}$ when their dimensions are clear. For positive definite (semi-
definite) matrices, we use the notation $P\succ 0$ $(P\succeq 0)$. For
negative definite (semi-definite) matrices, we use the notation $P\prec 0$
$(P\preceq 0)$. The notation $\text{col}[x_{1},\ldots,x_{n}]$ stands for the
column vector composed of the elements $x_{1},\ldots,x_{n}$. This notation is
also used when the components $x_{i}$ are vectors. The $\ell_{2}$ vector norm
(Euclidean norm) and the matrix norm induced by the $\ell_{2}$ vector norm are
both denoted as $||\cdot||$. We use $\mathcal{L}_{2}(0,T)$ (or simply
$\mathcal{L}_{2}$) to denote vector-valued functions
$z:[0,T]\to\mathbb{R}^{k}$ satisfying $\int_{0}^{T}\|z(t)\|^{2}dt<\infty$. For
a vector-valued signal $f(t)$ defined for all $t\geq 0$,
$||f(t)||_{\infty}:=\sup_{t\geq 0}||f(t)||$. For a differentiable function
$W:\mathbb{R}^{n}\to\mathbb{R}$ we denote by $\frac{\partial W}{\partial e}$
the row-vector of partial derivatives and by $\dot{W}(e)$ the total derivative
of $W(e)$ with respect to time.

##  II Problem Formulation

Consider the nonlinear system

$\left\\{\begin{aligned} \dot{x}(t)=&Ax(t)+Bu(t)+Sg(Vx(t),u(t))+F_{x}f(t),\\\
y(t)=&Cx(t)+F_{y}f(t),\end{aligned}\right.$ (1)

where $t\in\mathbb{R}^{+}$, $x\in{\mathbb{R}^{n}}$, $y\in{\mathbb{R}^{{m}}}$,
and $u\in{\mathbb{R}^{{l}}}$ are time, state, measured output and known input
vectors, respectively, ${n},{m},l\in\mathbb{N}$, matrices
$(A,B,S,V,C,F_{x},F_{y})$ of appropriate dimensions, and
$g:\mathbb{R}^{n_{v}}\times\mathbb{R}^{l}\to\mathbb{R}^{n_{g}}$ is a nonlinear
vector field. Function $f:\mathbb{R}^{+}\to\mathbb{R}^{n_{f}}$ denotes the
unknown fault vector, which contains both process and sensor faults,
$f_{x}(t)$ and $f_{y}(t)$, respectively,
$f(t):=\left[\begin{array}[]{ll}f_{x}^{T}(t)&f_{y}^{T}(t)\end{array}\right]^{T}$.
Matrix $F_{x}$ denotes the process fault distribution matrix while matrix
$F_{y}$ represents the contribution of the fault signal to sensor
measurements. Matrix $S$ indicates in which equation(s) the nonlinearity
appears explicitly, and $V$ is used to indicate which states play a role in
the nonlinearity. We often omit implicit time dependencies for notation
simplicity.

######  Assumption 1.

_(Globally Lipschitz Nonlinearity)_ Vector field $g(Vx,u)$ in (1) is globally
Lipschitz uniformly in $u(t)$ and $t$, i.e., there exists a known positive
constant $\alpha$ satisfying

$\|g(V\hat{x},u,t)-g(V{x},u,t)\|\leq\alpha\|V(\hat{x}-x)\|,$ (2)

for all $x,\hat{x}\in{\mathbb{R}^{n}}$, $u\in{\mathbb{R}^{l}}$, and
$t\in{\mathbb{R}^{+}}$.

######  Assumption 2.

_( $\mathcal{C}^{r}$ Fault Vector)_ The fault vector $f(t)$ in (1) is $r$
times differentiable, i.e., the time derivatives ${f}^{(1)}(t)$,
${f}^{(2)}(t)$, … ,$f^{(r)}(t)$ exist and are continuous, and $f^{(r)}(t)$ is
uniformly bounded.

###  II-A Ultra Local Fault Model

Under Assumption 2, we can write the linear approximated model
$\zeta_{i}(t)=f^{(i-1)}(t)$, $i=1,2,...,r-1$, $\zeta_{1}^{(r)}(t)=\mathbf{0}$,
for the fault vector as

$\left\\{\begin{aligned} \dot{\zeta}_{1}&=\zeta_{2},\\\
\dot{\zeta}_{2}&=\zeta_{3},\\\ \vdots\\\ \dot{\zeta}_{r-1}&=\zeta_{r},\\\
\dot{\zeta}_{r}&=\mathbf{0},\\\ f&=\zeta_{1}.\end{aligned}\right.$ (3)

Note that this model corresponds to an entry-wise $r$-th order Taylor time-
polynomial approximation at time $t$ of $f(t)$. The accuracy of the
approximated model increases as $f^{(r)}(t)$ goes to zero (entry-wise) and it
is exact for $f^{(r)}(t)=\mathbf{0}$. Model (3) is used to construct an
observer that ultra-locally [16]-[17] acts as a self-updating polynomial
spline approximating the actual value of the fault. To design such an
observer, we extend the system state, $x(t)$, with the states of the fault
model, $\zeta_{i}(t)$, $i=1,2,...,r$, and augment the system dynamics in (1)
with (3). We then design a Nonlinear Unknown Input Observer (NUIO) for the
augmented system to simultaneously estimate $x(t)$ and $f(t)$. We remark that
the number of the faults derivatives, $r$, added to the approximated model (3)
is problem-dependent and an optimal selection would depend on the frequency
characteristics of the fault. Increasing $r$ results in higher dimensional
augmented dynamics, and thus high dimensional observers as well. However,
having larger observers also provides more degrees of freedom for optimal
synthesis.

###  II-B Augmented Dynamics

Define the augmented state $x_{a}:=\text{col}[x,f,\dot{f},\ldots,f^{(r-1)}]$
and write the augmented dynamics using (1) and (3) as

$\displaystyle\left\\{\begin{aligned}
\dot{x}_{a}&=A_{a}x_{a}+B_{a}u+S_{a}g(V_{a}x_{a},u)+D_{a}f^{(r)},\\\
y&=C_{a}x_{a},\end{aligned}\right.$ (4a) where $\displaystyle A_{a}$
$\displaystyle:=\left[\begin{array}[]{ccccc}A&F_{x}&\boldsymbol{0}&\ldots&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{0}&I_{n_{f}}&\ldots&\boldsymbol{0}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\ldots&I_{n_{f}}\\\
\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\ldots&\boldsymbol{0}\end{array}\right],B_{a}:=\left[\begin{array}[]{cc}B\\\
\boldsymbol{0}\end{array}\right],$ (4b) $\displaystyle S_{a}$
$\displaystyle:=\left[\begin{array}[]{cc}S^{T}&\boldsymbol{0}\end{array}\right]^{T},V_{a}:=[\begin{array}[]{ll}V&\boldsymbol{0}\end{array}],$
$\displaystyle D_{a}$
$\displaystyle:=\left[\begin{array}[]{cc}\boldsymbol{0}&I_{n_{f}}\end{array}\right]^{T},C_{a}:=\left[\begin{array}[]{lll}C&F_{y}&\boldsymbol{0}\end{array}\right].$

###  II-C Joint State-Fault Nonlinear Observer

We propose the following nonlinear unknown input observer to estimate the
augmented state ${x}_{a}$

$\displaystyle\dot{z}$
$\displaystyle=Nz+Gu+Ly+MS_{a}g(V_{a}\hat{x}_{a}+J(y-C_{a}\hat{x}_{a}),u),$
(5a) $\displaystyle\hat{x}_{a}$ $\displaystyle=z-Ey,$ with observer state
$z\in{\mathbb{R}^{n_{z}}}$ ($n_{z}=n+rn_{f}$), estimate of the extended estate
$\hat{x}_{a}$, and matrices $(N,G,L,M)$ defined as $\displaystyle N$
$\displaystyle:=MA_{a}-KC_{a},$ $\displaystyle G$ $\displaystyle:=MB_{a},$
(5b) $\displaystyle L$ $\displaystyle:=K(I+C_{a}E)-MA_{a}E,$ $\displaystyle M$
$\displaystyle:=I+EC_{a}.$

Matrices $E$, $K$, and $J$ are observer gains to be designed. Note that the
fault estimate $\hat{f}(t)$ is given by $\hat{f}=\bar{C}\hat{x}_{a}$ with

$\bar{C}:={\left[\begin{array}[]{lll}\boldsymbol{0}_{n_{f}\times
n}&I_{n_{f}}&\boldsymbol{0}_{n_{f}\times n_{f}(r-1)}\end{array}\right]}.$ (6)

Define the estimation error as

$e:=\hat{x}_{a}-x_{a}=z-x_{a}-Ey=z-Mx_{a}.$

Then, the estimation error dynamics is given by

$\displaystyle\dot{e}$ $\displaystyle=Ne+(NM+LC_{a}-MA_{a})x_{a}+(G-MB_{a})u$
(7)
$\displaystyle+MS_{a}\big{(}g(V_{a}\hat{x}_{a}+J(y-C_{a}\hat{x}_{a}),u)-g(V_{a}{x}_{a},u)\big{)}$
$\displaystyle-MD_{a}f^{(r)}.$

Given the algebraic relations in (5b), it can be verify that
$G-MB_{a}=\mathbf{0}$ and $NM+LC_{a}-MA_{a}=\mathbf{0}$. Therefore, (7) can be
written as

$\displaystyle\dot{e}$ $\displaystyle=Ne+MS_{a}\delta g-MD_{a}f^{(r)},$ (8)

where $\delta
g:=g(V_{a}\hat{x}_{a}+J(y-C_{a}\hat{x}_{a}),u)-g(V_{a}{x}_{a},u)$. We have now
all the machinery required to state the problem we seek to solve.

######  Problem 1.

_(Fault Reconstruction)_ Consider the nonlinear system (1) with known input
and output signals, $u(t)$ and $y(t)$, and let Assumption 1 and Assumption 2
be satisfied. Further consider the approximated internal fault model (3), the
augmented dynamics (4a), and the observer (5). Design the observer gain
matrices $(E,K,J)$ so that: _1)_ all trajectories of the estimation error
dynamics (8) exists and are globally ultimately bounded uniformly in $t\geq
0$_;_ _2)_ $\int_{0}^{T}||\bar{C}e(t)||^{2}dt\leq
c\int_{0}^{T}\|f^{(r)}(t)\|^{2}dt$, for some known $c>0$, all
$T\in{\mathbb{R}^{+}}$, and $\bar{C}$ in (6); and _3)_ for
$f^{(r)}(t)=\mathbf{0}$, $t\geq 0$, $\lim_{t\rightarrow\infty}||e(t)||=0$.

Under Assumption 1 and Assumption 2, Problem 1 amounts to finding a fault
estimator that guarantees a bounded estimation error, $e(t)$, that $e(t)$ goes
to zero asymptotically if $f^{(r)}(t)=\mathbf{0}$, and that the
$\mathcal{L}_{2}(0,T)$ norm of $\bar{C}e(t)$ (the fault estimation error) is
upper bounded by that of $cf^{(r)}(t)$, for some $c>0$.

##  III Solution to Problem 1

###  III-A ISS Estimation Error Dynamics

In this section, we derive LMI conditions for designing the matrices $E,K$ and
$J$ of the augmented system observer in (5). As a stepping stone, we present a
sufficient condition for asymptotic stability of the origin of the estimation
error dynamics when there is no fault, then we prove boundedness of the
estimation error in the presence of faults using input-to-state stability
ideas [18].

######  Definition 1.

_(Input-to-State Stability[18, Def. 4.7])_ The error dynamics (8) is said to
be Input-to-State Stable (ISS) if there exist a class $\mathcal{K}\mathcal{L}$
function $\beta(\cdot)$ and a class $\mathcal{K}$ function $\gamma(\cdot)$
such that for any initial estimation error $e(t_{0})$ and any bounded
$f^{(r)}$, the solution $e(t)$ of (8) exists for all $t\geq t_{0}$ and
satisfies

$\|e(t)\|\leq\beta\left(\left\|e\left(t_{0}\right)\right\|,t-t_{0}\right)+\gamma\big{(}\|f^{(r)}(t)\|_{\infty}\big{)}.$
(9)

Note that ISS of the error dynamics (8) implies boundedness of the estimation
error for bounded $f^{(r)}(t)$. This follow directly from $\eqref{ISS_def}$.

######  Lemma 1.

_(ISS Lyapunov Function[18, Thm. 4.19])_ Consider the error dynamics (8) and
let $W(e)$ be a continuously differentiable function such that

$\alpha_{1}(\|e\|)\leq W(e)\leq\alpha_{2}(\|e\|),$ $\dot{W}(e)\leq-
W_{3}(e),\quad\forall\hskip 2.84526pt\|e\|\geq\xi(\|f^{(r)}\|),$

where $\alpha_{1}(\cdot)$ and $\alpha_{2}(\cdot)$ are class
$\mathcal{K}_{\infty}$ functions, $\xi(\cdot)$ is a class $\mathcal{K}$
function, and $W_{3}$ is a continuous positive definite function. Then, the
estimation error dynamics (8) is ISS with gain
$\gamma=\alpha_{1}^{-1}(\alpha_{2}(\xi))$.

Let $W(e):={e}^{T}P{e}$ be an ISS Lyapunov function candidate. Then, it
follows from (8) and the Lipschitz condition in (2) (i.e., $\|\delta
g\|\leq\alpha\|(V_{a}-JC_{a}){e}\|$) that

$\displaystyle\dot{W}(e)$ $\displaystyle\leq e^{T}\Delta
e-2e^{T}PMD_{a}f^{(r)},$ (10)

where $\Delta:=N^{T}P+PN+\alpha
PMS_{a}S_{a}^{T}M^{T}P+\alpha(V_{a}-JC_{a})^{T}(V_{a}-JC_{a})$. A complete
derivation of inequality (10) is given in Appendix A-A. Inequality (10)
implies the following

$\displaystyle\dot{W}(e)\leq$
$\displaystyle-\lambda_{\min}(-\Delta)\|e\|^{2}+2\|e\|\|PMD_{a}\|\|f^{(r)}\|$
(11) $\displaystyle=$
$\displaystyle-(1-\theta)\lambda_{\min}(-\Delta)\|e\|^{2}-\theta\lambda_{\min}(-\Delta)\|e\|^{2}$
$\displaystyle+2\|e\|\|PMD_{a}\|\|f^{(r)}\|,$

for any $\theta\in(0,1)$. Therefore, by (11) and Lemma 1, if $\Delta$ is
negative definite, system (8) is ISS with input $f^{(r)}$ and linear ISS-gain

$\gamma(\|f^{(r)}\|)=\frac{2\|PMD_{a}\|}{\theta\lambda_{\min}(-\Delta)}\|f^{(r)}\|.$
(12)

Based on the above discussion, the next proposition formalizes a LMI condition
($\Delta\prec 0$) that guarantees an ISS estimation error dynamics (8).
Without loss of generality, for numerical tractability, we enforce
$\Delta+\epsilon I\preceq 0$ for some small given $\epsilon>0$ instead of
$\Delta\prec 0$.

######  Proposition 1.

_(ISS Estimation Error Dynamics)_ Consider the error dynamics (8) and suppose
there exist matrices ${\mathbb{R}^{n_{z}\times n_{z}}}\ni P\succ
0,R\in{\mathbb{R}^{n_{z}\times m}},Q\in{\mathbb{R}^{n_{z}\times m}}$, and
$J\in{\mathbb{R}^{n_{v}\times m}}$ satisfying the inequality

$\left[\begin{array}[]{cc}X+\epsilon I&X_{12}\\\
X_{12}^{T}&-I\end{array}\right]\preceq 0,$ (13)

for some given $\epsilon>0$ and matrices $X$ and $X_{12}$ defined as

$\displaystyle X:=$ $\displaystyle
S_{11}+\alpha(V_{a}^{T}V_{a}-V_{a}^{T}JC_{a}-C_{a}^{T}J^{T}V_{a}),$ (14)
$\displaystyle S_{11}:=$ $\displaystyle
A_{a}^{\top}P+A_{a}^{\top}C_{a}^{\top}R^{\top}-C_{a}^{T}Q^{T}+PA_{a}$
$\displaystyle+RC_{a}A_{a}-QC_{a},$ $\displaystyle X_{12}:=$
$\displaystyle\sqrt{\alpha}{\left[\begin{array}[]{ll}(P+RC_{a})S_{a}&C_{a}^{T}J^{T}\end{array}\right],}$

with $\alpha$ from (2) and the remaining matrices in (4b); then, the ISS-gain
from input $f^{(r)}$ to the estimation error $e$ in (8) is upper bounded by
${2\|P(I+EC_{a})D_{a}\|}\epsilon^{-1}$ with $E=P^{{}^{-1}}R$.

_Proof_ : We want to prove that $\Delta+\epsilon I\preceq 0$ is equivalent to
(13)-(14). Consider the expression for $\Delta$ in the text below (10). This
expression is written in terms of the original observer gains $(E,K,J)$. Using
this $\Delta$, (5b), and Schur complements on $\Delta+\epsilon I\preceq 0$, we
can write (13) with $X$ and $X_{12}$ in terms of the original observer gains
$(E,K,J)$ as

$\displaystyle X$
$\displaystyle:=N^{T}P+PN+\alpha(V_{a}^{T}V_{a}-V_{a}^{T}JC_{a}-C_{a}^{T}J^{T}V_{a}),$
$\displaystyle X_{12}$
$\displaystyle:=\sqrt{\alpha}{\left[\begin{array}[]{ll}PMS_{a}&C_{a}^{T}J^{T}\end{array}\right]},$

where $N^{T}P+PN$ expands as

$A_{a}^{T}\left(I+EC_{a}\right)^{T}P-C_{a}^{T}KP+P\left(I+EC_{a}\right)A_{a}-PKC_{a},$

and $PMS_{a}$ becomes

$PS_{a}+PEC_{a}S_{a}.$

Consider the following change of variables

$R:=PE,\qquad Q:=PK.$ (15)

Applying (15) on the above expanded $X$ and $X_{12}$, the linear inequality
(13)-(14) can be concluded. Clearly, $\Delta+\epsilon I\preceq 0$ implies
$\lambda_{\min}(-\Delta)\geq\epsilon$. Then, using (12), we can conclude the
bound on the ISS-gain. $\blacksquare$

######  Remark 1.

_(LMI Feasibility)_ If $F_{y}$ is full row rank (i.e., there as many sensor
faults, $f_{y}$, as sensors), the LMI in (13) is always infeasible
(observability is lost) [19, 20, 21]. The standard practice to circumvent this
issue is to assume $\text{rank}[F_{y}]<m$ [19, 20, 21]. Note that this is only
a necessary condition for the estimator to exist, but it does not guarantee
the LMI in (13) to be feasible (this has to be checked on a case by case
basis).

######  Remark 2.

_(Decoupling of Nonlinearities)_ Note that it might be possible to cancel the
effect of the nonlinearity $\delta g$ in the estimation error dynamics (8)
provided that we can select $M$ to satisfy $MS_{a}=\mathbf{0}$. This algebraic
condition can be written in terms of the observer gain $E$ using the
definition of $M$ in (5b) as follows

$EC_{a}S_{a}=-S_{a}.$ (16)

Equation (16) has a solution $E$, if and only if matrix $C_{a}S_{a}$ has full
column rank (details about this can be found in [22]). Physically, the latter
means there must be measurements from which the nonlinearity appears in the
output dynamics. For the sake of generality, we do not assume this thought. We
tackle the problem without imposing $MS_{a}=\mathbf{0}$ (so considering the
nonlinear terms in the estimation error dynamics).

Proposition 1 is used to enforce that all trajectories of the estimation error
dynamics (8) exist and are bounded for all $t\geq 0$. If the observer gains
$(E,K,J)$ satisfy the ISS LMI in (13), boundedness in guaranteed. This follows
directly from Assumption 2 and Definition 1.

###  III-B $\mathcal{L}_{2}$ Performance Criteria

To maximize the performance of the reconstruction scheme, we seek to minimize
the effect of $f^{(r)}$ (treated as an external disturbance) on the estimation
error dynamics (8). We could use the ISS formulation in Proposition 1 to cast
an optimization problem where we minimize the ISS gain and treat the LMI in
(13) as an optimization constraint. By doing so, we would be reducing the
effect of $f^{(r)}$ on the complete vector of estimation errors $e$ (state,
fault, and fault derivatives estimation errors). Note, however, that the
purpose of the filter is to reconstruct faults only, so the performance in
state estimation and in the error of higher order fault derivatives is not
relevant. To this end, we seek to minimize the $\mathcal{L}_{2}$ gain from
$f^{(r)}$ to the fault estimation error $\bar{C}e(t)$, with $\bar{C}$ in (6).
The $\mathcal{L}_{2}$ gain allows for an input-output
($f^{(r)}\rightarrow\bar{C}e(t)$) characterization of performance.

######  Definition 2.

_( $\mathcal{L}_{2}$-gain [23])_ We say that the estimation error dynamics (8)
with input $f^{(r)}$ and output $\bar{C}e(t)$ has a $\mathcal{L}_{2}$ gain
less than or equal to $\lambda$ if the following is satisfied

$\int_{0}^{T}\|\bar{C}e(t)\|^{2}dt\leq\lambda^{2}\int_{0}^{T}\|f^{(r)}(t)\|^{2}dt,$

for all $T\geq 0$ and $f^{(r)}(t)\in L_{2}(0,T)$.

In the following lemma, we state a Lyapunov-based sufficient condition (the
Hamilton-Jacobi inequality) for having a bounded $\mathcal{L}_{2}$ gain (see
[18, Thm. 5.5] and [23, Thm. 2] for further details).

######  Lemma 2.

_( $\mathcal{L}_{2}$-gain Inequality [18, Thm. 5.5])_ Consider (8) and suppose
there exists a continuously differentiable positive semi-definite function
$W(e)$ satisfying

$\displaystyle\frac{\partial W}{\partial e}(Ne+MS_{a}\delta
g)+\frac{1}{2\lambda^{2}}\frac{\partial W}{\partial
e}MD_{a}D_{a}^{T}M^{T}\left(\frac{\partial W}{\partial e}\right)^{T}$ (17)
$\displaystyle\hskip 14.22636pt+\frac{1}{2}p^{T}p\leq 0,$

with $p=\bar{C}e$ and $\bar{C}$ as in (6). Then, the $\mathcal{L}_{2}$-gain
from $f^{(r)}$ to the fault estimation error $p=\bar{C}e$ in (8) is less than
or equal to $\lambda$.

Based on the Lemma 2, the next proposition formalizes an LMI-based condition
guaranteeing the inequality (17) to have the finite $\mathcal{L}_{2}$-gain of
the mapping from $f^{(r)}$ to the fault estimation error.

######  Proposition 2 (Finite $\mathcal{L}_{2}$-gain).

Suppose there exist matrices $P\succ 0,R,$ $Q$, $J$ and scalar $\rho\geq 0$
satisfying

$\displaystyle\left[\begin{array}[]{ccc}L_{11}&(P+RC_{a})D_{a}&X_{12}\\\
&-{\rho}I&\boldsymbol{0}\\\ &*&-I\end{array}\right]\preceq 0,$ (21)
$\displaystyle L_{11}:=X+\frac{1}{2}\bar{C}^{T}\bar{C},$ (22)

and $X,X_{12}$ as defined in (14), $\bar{C}$ in (6) and the remaining matrices
in (4b). Then, the $\mathcal{L}_{2}$-gain from $f^{(r)}$ to the fault
estimation error $p=\bar{C}e$ in (8) is upper bounded by $\sqrt{2\rho}$.

_Proof :_ The proof is similar to the proof of Proposition 1 and is given in
Appendix A-B. $\blacksquare$

Using Proposition 2, we next cast a semidefinite program where we seek to
minimize the $\mathcal{L}_{2}$-gain from $f^{(r)}$ to $p=\bar{C}e$. We add the
ISS LMI in (13) as a constraint to this program to enforce that the resulting
filter also guarantees boundedness for bounded faults. The latter is important
to avoid that the filter diverges (as the $\mathcal{L}_{2}$ criteria does not
guarantee stability).

######  Theorem 1.

_(Optimal Fault Estimator)_ Consider the augmented system dynamics (4), the
filter (5), and the corresponding estimation error dynamics (8). To design the
optimal fault estimator (5), solve the following convex program

$\begin{array}[]{cl}\min\limits_{P,R,Q,J,\rho}&\rho\\\ \text{s.t.}&\vspace{1
mm}\left[\begin{array}[]{cc}X+\epsilon I&X_{12}\\\
X_{12}^{T}&-I\end{array}\right]\preceq 0\\\
&\left[\begin{array}[]{ccc}L_{11}&(P+RC_{a})D_{a}&X_{12}\\\
&-{\rho}I&\boldsymbol{0}\\\ &*&-I\end{array}\right]\preceq 0\\\\[14.22636pt]
&P\succ 0,\qquad\rho\geq 0,\end{array}$ (23)

with given $\epsilon>0$, $X,X_{12}$ as defined in (14), $L_{11}$ in (22), and
the remaining matrices in (4b). Let us denote the optimizers by
$P^{\star},R^{\star},Q^{\star},J^{\star},\rho^{\star},$ and based on (15),
define the matrices
$E^{\star}=P^{\star^{-1}}R^{\star},K^{\star}=P^{\star^{-1}}Q^{\star}$. Then,
we have the following:

  1. 1.

The ISS-gain from $f^{(r)}$ to $e$ is upper bounded by
${2\|P^{\star}(I+E^{\star}C_{a})D_{a}\|}\epsilon^{-1}$.

  2. 2.

The $\mathcal{L}_{2}$-gain from $f^{(r)}$ to $\bar{C}e$ is upper bounded by
$\sqrt{2\rho^{\star}}$.

_Proof :_ Theorem 1 follows from the above discussion, Proposition 1, and
Proposition 2. $\blacksquare$

######  Remark 3 (Perfect Estimation).

The developed methodology can guarantee zero estimation error for zero
$f^{(r)}$, i.e., when the $r$-time derivative of the fault vector vanishes.
This follows directly from the estimation error dynamics (8) since if
$f^{(r)}=\boldsymbol{0}$, the error dynamics becomes independent of external
disturbances.

Figure 1: Benchmark System Schematic.

##  IV Simulation Results

In this section, we evaluate the proposed method by a benchmark example for
FDI [24, 25, 26]. The dynamics of the system (a single-link robotic arm with a
revolute elastic joint, see Fig. 1 for a schematic) can be described as
follows:

$\displaystyle
J_{l}\ddot{q}_{l}+F_{l}\dot{q}_{l}+k\left(q_{l}-q_{m}\right)+mgc\sin\left(q_{l}\right)=0,$
$\displaystyle
J_{m}\ddot{q}_{m}+F_{m}\dot{q}_{m}-k\left(q_{l}-q_{m}\right)=k_{\tau}u,$

where ${q}_{l}$ and ${q}_{m}$ are the angular position of the link and the
angular position of the motor, respectively. Constants $J_{l}$ and $J_{m}$ are
the moments of inertia of the link and the motor and $F_{l}$ and $F_{m}$ are
the viscous coefficients associated to friction acting at the link and the
motor, respectively. The flexibility in the joint is modeled by a spring with
a spring coefficient $k$, $m$ is the link mass, $g$ is the gravity constant,
$c$ is the height of the link center of mass, $k_{\tau}$ is the amplifier
gain, and $u$ is the torque input delivered by the motor. Units are in SI, and
the parameters are: $J_{l}=4.5,J_{m}=1,F_{l}=0.5,F_{m}=1,k=2,m=4,g=9.8,c=0.5,$
and $k_{\tau}=1$. The torque input is $u=2\text{sin}(0.25t)$.

By selecting $x_{1}:=\dot{q}_{m},x_{2}:=q_{m},x_{3}:=\dot{q}_{l}$, and
$x_{4}:=q_{l}$, the system can be written in the form of (1):

$\displaystyle\dot{x}$ $\displaystyle=Ax+Bu+Sg\left(Vx\right)+F_{x}f$ (24)
$\displaystyle y$ $\displaystyle=Cx+F_{y}f,$

where $x:=[x_{1},x_{2},x_{3},x_{4}]^{T}$ is the state vector, and

$\displaystyle A$
$\displaystyle=\left[\begin{array}[]{cccc}-\frac{F_{m}}{Jm}&-\frac{k}{J_{m}}&0&\frac{k}{J_{m}}\\\
1&0&0&0\\\ 0&\frac{k}{J_{l}}&-\frac{F_{l}}{J_{l}}&-\frac{k}{J_{l}}\\\
0&0&1&0\end{array}\right],$ $\displaystyle B$
$\displaystyle=\left[\begin{array}[]{cccc}\frac{k_{\tau}}{J_{m}}&0&0&0\\\
\end{array}\right]^{T},$ $\displaystyle C$
$\displaystyle=\left[\begin{array}[]{cccc}0&1&0&0\\\
0&0&0&1\end{array}\right],$ $\displaystyle S$
$\displaystyle=\left[\begin{array}[]{cccc}0&0&\frac{mgc}{J_{l}}&0\\\
\end{array}\right]^{T},$ $\displaystyle V$
$\displaystyle=\left[\begin{array}[]{cccc}0&0&0&1\\\ \end{array}\right].$

The nonlinearity is given by $g\left(Vx\right)=-\sin\left(x_{4}\right)$, which
is Lipschitz with constant $\alpha=1$. We set initial conditions as
$x(0)=[0.01,0.01,0.01,0.01]^{T}$. In addition, the first measurement is the
angular position of the motor, and the second one is the angular position of
the link. Then, we simulate for two different scenarios, one for the first
sensor fault where we have

$\displaystyle F_{x}$ $\displaystyle=\left[\begin{array}[]{c}0\\\ 0\\\ 0\\\
0\\\ \end{array}\right],$ $\displaystyle F_{y}$
$\displaystyle=\left[\begin{array}[]{cc}1\\\ 0\\\ \end{array}\right],$

and the other one for the actuator fault with

$\displaystyle F_{x}$ $\displaystyle=\left[\begin{array}[]{c}1\\\ 0\\\ 0\\\
0\\\ \end{array}\right],$ $\displaystyle F_{y}$
$\displaystyle=\left[\begin{array}[]{cc}0\\\ 0\\\ \end{array}\right].$

For both scenarios, we augment the benchmark system (24) using (4) with $r=1$.
Next, we design an NUIO observer of the form (5) by solving the minimization
problem (23) in Theorem 1, for the estimation of the actuator fault and the
first sensor fault. To solve the convex program (23), we use the YALMIP
toolbox in MATLAB.

Figure 2: The actual sensor fault and its estimate using the proposed method.

Finally, we apply the designed NUIO to the system in the simulation. The
initial condition of the observer in the simulation is taken as the zero
vector. A gradually increasing (similar to an incipient fault [27]) sensor
fault occurs at time $t=25$ $(sec)$. Figure 2 depicts the estimated actuator
fault and its actual value. For actuator fault, to indicate the capability of
the developed method to estimate time-varying faults, a sinusoidal actuator
fault with the same frequency of input is simulated (i.e.,
$f_{x}=0.1sin(0.25t)$). Figure 3 shows the estimated actuator fault and its
actual value. It can be seen that the estimated actuator and first sensor
faults follow the actual fault with almost zero error in steady state.

##  V Conclusion

In this paper, a method for time-varying actuator and sensor faults estimation
in nonlinear systems has been proposed. To this end, we have augmented the
system dynamics with an approximated internal model for the fault vector
(using ideas from ultra-local models). Then, a nonlinear unknown input
observer is constructed for the augmented dynamics. The fault estimation
problem has been re-formulated as a robust state estimation problem in the
augmented state. We have provide sufficient conditions that ensure the
stability of the observer error dynamics (guaranteeing asymptotic fault
estimation) and robustness against model mismatch in the internal fault model
(in terms of a finite $\mathcal{L}_{2}$-gain from fault model mismatch to
fault estimation error). In Theorem 1, a convex minimization problem has been
casted for destining optimal estimators. Finally, simulations for a benchmark
system are given to illustrate the performance of the proposed methodology.
The numerical simulations show that the proposed method properly estimates
sensor and actuator faults for the benchmark system.

Figure 3: The actual actuator fault and its estimate using the proposed
method.

##  Appendix A Appendices

###  A-A Lyapanov Function Proof

The inequality in (10) can be derived as follows:

$\displaystyle\dot{W}=$ $\displaystyle
e^{T}\left(N^{T}P+PN\right)e+2e^{T}PMS_{a}\delta g$ (25)
$\displaystyle-2e^{T}PMD_{a}f^{(r)}$ $\displaystyle\leq$ $\displaystyle
e^{T}\left(N^{T}P+PN\right)e+2\left\|e^{T}PMS_{a}\right\|\|\delta g\|$
$\displaystyle-2e^{T}PMD_{a}f^{(r)}$ $\displaystyle\leq$ $\displaystyle
e^{T}\left(N^{T}P+PN\right)e$
$\displaystyle+2\left\|e^{T}PMS_{a}\right\|\alpha\|(V_{a}-JC_{a}){e}\|-2e^{T}PMD_{a}f^{(r)}$
$\displaystyle\leq$ $\displaystyle
e^{T}\left(N^{T}P+PN\right)e+\alpha\big{(}\left\|e^{T}PMS_{a}\right\|^{2}+$
$\displaystyle\|(V_{a}-JC_{a})e\|^{2}\big{)}-2e^{T}PMD_{a}f^{(r)}$
$\displaystyle=$ $\displaystyle e^{T}\big{(}N^{T}P+PN+\alpha
PMS_{a}S_{a}^{T}M^{T}P$
$\displaystyle+\alpha(V_{a}-JC_{a})^{T}(V_{a}-JC_{a})\big{)}e-2e^{T}PMD_{a}f^{(r)}$
$\displaystyle:=$ $\displaystyle e^{T}\Delta e-2e^{T}PMD_{a}f^{(r)}.$

###  A-B $\mathcal{L}_{2}$-Gain LMI

If we define $W(e):={e}^{T}P{e}$ with positive definite matrix $P$; then (17)
can be written as follows

$\displaystyle
e^{T}\left(N^{T}P+PN+\frac{1}{2}\bar{C}^{T}\bar{C}\right)e+2e^{T}PMS_{a}\delta
g$ $\displaystyle+\frac{2}{\lambda^{2}}e^{T}PMD_{a}D_{a}^{T}M^{T}Pe\leq 0.$

We can find an upper bound for the above inequality, similar to stability
proof in Appendix A-A, as

$e^{T}(\Delta+\frac{1}{2}\bar{C}^{T}\bar{C}+\frac{2}{\lambda^{2}}PMD_{a}D_{a}^{T}M^{T}P)e\leq
0,$

where $\Delta$ is same as definition in (25). Using Schur complement the above
inequality is equivalent to

$\left[\begin{array}[]{cc}\Delta+\frac{1}{2}\bar{C}^{T}\bar{C}&PMD_{a}\\\
&-\frac{\lambda^{2}}{2}I\end{array}\right]\preceq 0.$

If we follow the same procedure in proof of Proposition 1, the equivalent
inequality can be found as

$\left[\begin{array}[]{ccc}X+\frac{1}{2}\bar{C}^{T}\bar{C}&(P+RC_{a})D_{a}&X_{12}\\\
&-\frac{\lambda^{2}}{2}I&0\\\ &*&-I\end{array}\right]\preceq 0.$

Finally, by defining a change of variable as $\rho:=\frac{\lambda^{2}}{2}$,
(21) can be implied.

## acknowledgment

This publication is part of the project Digital Twin project 4.3 with project
number P18-03 of the research programme Perspectief which is (mainly) financed
by the Dutch Research Council (NWO).

## References

  * [1] S. X. Ding, _Model-based fault diagnosis techniques: design schemes, algorithms, and tools_. Springer Science & Business Media, 2008.
  * [2] X. Zhang, M. M. Polycarpou, and T. Parisini, “A robust detection and isolation scheme for abrupt and incipient faults in nonlinear systems,” _IEEE Transactions on Automatic Control_ , vol. 47, no. 4, pp. 576–593, 2002.
  * [3] I. Hwang, S. Kim, Y. Kim, and C. E. Seah, “A survey of fault detection, isolation, and reconfiguration methods,” _IEEE Transactions on Control Systems Technology_ , vol. 18, no. 3, pp. 636–653, 2009.
  * [4] A. Moradmand, B. Shafai, and M. Saif, “A design procedure for robust actuator and sensor fault detection,” in _2020 7th International Conference on Control, Decision and Information Technologies (CoDIT)_ , vol. 1. IEEE, 2020, pp. 709–714.
  * [5] J. Chen and R. J. Patton, _Robust model-based fault diagnosis for dynamic systems_. Springer Science & Business Media, 2012, vol. 3.
  * [6] S. Narasimhan, P. Vachhani, and R. Rengaswamy, “New nonlinear residual feedback observer for fault diagnosis in nonlinear systems,” _Automatica_ , vol. 44, no. 9, pp. 2222–2229, 2008.
  * [7] P. Mohajerin Esfahani and J. Lygeros, “A tractable fault detection and isolation approach for nonlinear systems with probabilistic performance,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 3, pp. 633–647, 2015\.
  * [8] C. De Persis and A. Isidori, “A geometric approach to nonlinear fault detection and isolation,” _IEEE Transactions on Automatic Control_ , vol. 46, no. 6, pp. 853–865, 2001.
  * [9] M. Liu and P. Shi, “Sensor fault estimation and tolerant control for itô stochastic systems with a descriptor sliding mode approach,” _Automatica_ , vol. 49, no. 5, pp. 1242–1250, 2013.
  * [10] M. Liu, X. Cao, and P. Shi, “Fuzzy-model-based fault-tolerant design for nonlinear stochastic systems against simultaneous sensor and actuator faults,” _IEEE Transactions on Fuzzy Systems_ , vol. 21, no. 5, pp. 789–799, 2012.
  * [11] B. Jiang, M. Staroswiecki, and V. Cocquempot, “Fault accommodation for nonlinear dynamic systems,” _IEEE Transactions on Automatic Control_ , vol. 51, pp. 1578–1583, 2006.
  * [12] C. P. Vo, H. V. Dao, K. K. Ahn, _et al._ , “Robust fault-tolerant control of an electro-hydraulic actuator with a novel nonlinear unknown input observer,” _IEEE Access_ , vol. 9, pp. 30 750–30 760, 2021.
  * [13] J.-W. Zhu, G.-H. Yang, H. Wang, and F. Wang, “Fault estimation for a class of nonlinear systems based on intermediate estimator,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 9, pp. 2518–2524, 2015.
  * [14] M. HOU and R. PATTON, “Input observability and input reconstruction,” _Automatica_ , vol. 34, pp. 789–794, 1998.
  * [15] C. van der Ploeg, M. Alirezaei, N. van de Wouw, and P. Mohajerin Esfahani, “Multiple faults estimation in dynamical systems: Tractable design and performance bounds,” _arXiv preprint arXiv:2011.13730_ , 2020.
  * [16] M. Fliess and C. Join, “Model-free control,” _International Journal of Control_ , vol. 86, pp. 2228–2252, 2013.
  * [17] H. Sira-Ramírez, A. Luviano-Juárez, M. Ramírez-Neria, and E. W. Zurita-Bustamante, _Active disturbance rejection control of dynamic systems: a flatness based approach_. Butterworth-Heinemann, 2018.
  * [18] H. K. Khalil, “Nonlinear systems third edition,” _Patience Hall_ , vol. 115, 2002.
  * [19] C. Edwards and C. P. Tan, “A comparison of sliding mode and unknown input observers for fault reconstruction,” _European Journal of Control_ , vol. 12, no. 3, pp. 245–260, 2006.
  * [20] C. Edwards and S. K. Spurgeon, “A sliding mode observer based fdi scheme for the ship benchmark,” _European Journal of Control_ , vol. 6, no. 4, pp. 341–355, 2000.
  * [21] P. Frank and J. Wunnenberg, “Robust fault diagnosis using unknown input observer schemes,” _Fault Diagnosis in Dynamic Systems_ , pp. 47–97, 1989\.
  * [22] W. Chen and M. Saif, “Unknown input observer design for a class of nonlinear systems: an lmi approach,” in _American Control Conference_. IEEE, 2006, pp. 5–pp.
  * [23] A. J. Van Der Schaft, “L 2-gain analysis of nonlinear systems and nonlinear state feedback $h_{\infty}$ control,” _IEEE Transactions on Automatic Control_ , vol. 37, no. 6, pp. 770–784, 1992.
  * [24] C. Keliris, M. M. Polycarpou, and T. Parisini, “An integrated learning and filtering approach for fault diagnosis of a class of nonlinear dynamical systems,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 28, no. 4, pp. 988–1004, 2016.
  * [25] V. Reppa, M. M. Polycarpou, and C. G. Panayiotou, “Adaptive approximation for multiple sensor fault detection and isolation of nonlinear uncertain systems,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 25, no. 1, pp. 137–153, 2013.
  * [26] X. Zhang, T. Parisini, and M. M. Polycarpou, “Sensor bias fault isolation in a class of nonlinear systems,” _IEEE Transactions on Automatic Control_ , vol. 50, no. 3, pp. 370–376, 2005.
  * [27] X. Zhang, M. M. Polycarpou, and T. Parisini, “A robust detection and isolation scheme for abrupt and incipient faults in nonlinear systems,” _IEEE Transactions on Automatic Control_ , vol. 47, no. 4, pp. 576–593, 2002.


# Beware the recent past: a bias in spectral energy distribution modeling due
to bursty star formation

P. Haskell1, S. Das1, D. J. B. Smith${}^{1},$ R. K. Cochrane${}^{2,3},$ C. C.
Hayward3 and D. Anglés-Alcázar${}^{4}{}^{,}{}^{3}$
1Centre for Astrophysics Research, University of Hertfordshire, College Lane,
Hatfield AL10 9AB, UK
2Department of Astronomy, Columbia University, New York, NY 10027, USA
3Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue,
New York, NY 10010, USA
4Department of Physics, University of Connecticut, 196 Auditorium Road,
U-3046, Storrs, CT 06269-3046, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We investigate how the recovery of galaxy star formation rates (SFRs) using
energy-balance spectral energy distribution (SED) fitting codes depends on
their recent star formation histories (SFHs). We use the Magphys and
Prospector codes to fit 6,706 synthetic spectral energy distributions (SEDs)
of simulated massive galaxies at $1<z<8$ from the Feedback in Realistic
Environments (FIRE) project. We identify a previously unknown systematic error
in the Magphys results due to bursty star formation: the SFR estimates of
individual galaxies can differ from the true values by as much as 1 dex, at
large statistical significance ($>5\sigma$), depending on the details of their
recent SFH. The SFRs inferred using Prospector do not exhibit this trend,
likely because unlike Magphys, Prospector uses non-parametric SFHs. We urge
caution when using Magphys, or other codes assuming parametric SFHs, to study
galaxies where the average SFR may have changed significantly over the last
$\sim$100 Myr, such as those which have recently quenched their star formation
or those experiencing an ongoing burst. This concern is especially relevant,
for example, when fitting JWST observations of very high-redshift galaxies.

###### keywords:

galaxies: fundamental parameters - methods: data analysis - galaxies:
photometry

††pubyear: 2023††pagerange: Beware the recent past: a bias in spectral energy
distribution modeling due to bursty star formation–Beware the recent past: a
bias in spectral energy distribution modeling due to bursty star formation

##  1 Introduction

Star formation rates (SFRs) and star formation histories (SFHs) are both of
great importance to our understanding of galaxies. SFRs are critical since
(along with their stellar mass) they enable us to put galaxies in their
contemporary context relative to the evolving SFR–stellar mass relation (e.g.
Noeske et al., 2007). SFHs encode the build-up of stellar mass in the Universe
(Dye, 2008), potentially yielding insight into the physical processes that
drive galaxy evolution (e.g. Madau et al., 1996; Behroozi et al., 2013; Wang
et al., 2022). SFRs and other galaxy properties are typically inferred for
high-redshift galaxies by fitting a model spectral energy distribution (SED)
to the available data, using whatever photometry is available. There are many
freely-available SED fitting codes able to do this (as discussed by e.g.
Walcher et al., 2011; Conroy, 2013; Baes, 2019).

SED models frequently differ in their approach and assumptions in important
ways which potentially impact the fidelity of inferred galaxy properties (Hunt
et al. 2019, Pacifici et al. 2023). For example, they may employ different
stellar population synthesis models (e.g. Bruzual & Charlot 2003, Maraston
2005, Conroy et al. 2009; see also Baldwin et al. 2018), initial mass
functions (IMFs; e.g. Salpeter 1955, Kroupa 2001, Chabrier 2003; see also
Kroupa & Jerabkova 2019), dust models (e.g. Calzetti 1997, Charlot & Fall
2000, Draine & Li 2007, da Cunha et al. 2008) and SFH forms. These can be
parametric or non-parametric; while the former typically use one of many
simple functional forms to represent the SFH (e.g. Behroozi et al., 2013;
Simha et al., 2014; Carnall et al., 2019), the latter can be implemented e.g.
by partitioning the history into time bins and allocating stellar mass to each
bin from a chosen prior distribution as the observations allow (e.g. Leja et
al. 2019).

Many SED fitters, including two of the most well-known codes – Magphys (da
Cunha et al., 2008, 2015) and Prospector (Leja et al., 2017) – assume equality
(or balance) between the energy absorbed by dust from starlight at short
wavelengths and that which is re-radiated in the far-infrared. Physically,
this must be true for the total absorbed and re-radiated luminosities (i.e.
integrating over solid angle), but strict equality should not hold for every
line of sight due to the non-isotropic nature of dust attenuation. Recent
observations of spatial offsets between the starlight and dust in high-
redshift galaxies (e.g. Hodge et al., 2016; Rujopakarn et al., 2019; Cochrane
et al., 2021) cast doubt on whether these codes are appropriate in such cases
where there is the potential for a lack of energy balance along a particular
line of sight, to the point that other fitters explicitly avoid this
assumption (e.g. Liu, 2020; Kokorev et al., 2021). However, in Haskell et al.
(2023, hereafter H23), we showed that the performance of Magphys is similar
for galaxies with UV/FIR spatial offsets as it is for local galaxies (Hayward
& Smith, 2015), suggesting that energy-balance codes remain the gold standard
for deriving the properties of distant galaxies from photometric surveys.

With myriad codes available, with different physics and underlying
assumptions, it is critical to validate the performance of each method. As
well as the works of this group, which have compared the physical properties
inferred using energy-balance SED fitting with Magphys to simulated galaxies
for which the ‘ground truth’ properties are known (Hayward & Smith 2015; Smith
& Hayward 2015, 2018; H23), there is a large body of work by other authors
aiming to understand how well SED fitting works in different ways (e.g. da
Cunha et al. 2008, Lee et al. 2009, Noll et al. 2009, Wuyts et al. 2009, Lee
2010, Pforr et al. 2012, Simha et al. 2014, Mobasher et al. 2015, Hunt et al.
2019, Dudzevičiūtė et al. 2020, Lower et al. 2020, Pacifici et al. 2023, Best
et al. 2023). Since validation tests are often accomplished by averaging over
a large population of galaxies, there is always danger of “washing out”
important details or systematic effects that could impact our understanding of
how well SED fitting methods work in real-Universe situations.

In this Letter we build on the analysis of H23 and run the Prospector SED
fitting code on the same synthetic galaxy SEDs for which we already have
Magphys fits, focusing on how the fidelity of the inferred SFR depends on the
recent (i.e. within the last 100 Myr) SFH. Galaxies experiencing rapid changes
in SFR over short ($\sim$10s-100s of Myr) timescales have been identified both
observationally (e.g. Guo et al. 2016, Faisst et al. 2019, Broussard et al.
2022, Looser et al. 2023) and in simulations (e.g. Sparre et al. 2017, Iyer et
al. 2020, Flores Velázquez et al. 2021, Hopkins et al. 2023).

In the FIRE simulations, massive galaxies at high redshift and dwarfs at all
redshifts exhibit very bursty SFHs, with short ($\sim$10-Myr) bursts followed
by temporarily quenched periods111We will use the shorthand ‘temp-quenched’ to
refer to the periods of low-specific SFR in between bursts that are common in
the FIRE simulations analysed here. The term ‘mini-quenched’ has been used by
some authors, but given that ‘mini’ refers to size rather than temporal
duration, we believe that ‘temp-quenched’ is more accurate. We use the term
‘starbursting’ to refer to simulated galaxies that are currently undergoing a
burst, i.e. the 10 Myr-averaged SFR is greater than the 100 Myr-averaged
value. lasting a few 100 Myr; see Sparre et al. (2017) for a detailed analysis
and Hayward & Hopkins (2017) and Hopkins et al. (2023) for discussions of the
physical origin of this behavior. They thus provide an excellent test data set
for exploring the effects of bursty star formation on SED modeling
inferences.222Owing to the lack of AGN feedback in these simulations, none of
the simulated galaxies are permanently quenched. However, because the bias
depends on the shape of the SFH over the past 100 Myr, our conclusions should
apply equally well to galaxies that are undergoing quenching and will remain
quenched permanently, not only ‘temp-quenched’ galaxies.

We demonstrate that Magphys results exhibit a systematic bias in the inferred
SFR due to bursty star formation: the SFR is overestimated (underestimated)
for simulated galaxies experiencing a burst (temporarily quenched period). We
show that Prospector does not exhibit this bias, likely due to its use of non-
parametric SFHs.

The remainder of this Letter is structured as follows. In Section 2 we outline
the tools and methods used to create the observations and to fit the SEDs,
while in Section 3 we present the results of the fitting, and in Section 4
discuss the implications. In Section 5 we make some concluding remarks. We
adopt a standard cosmology with $H_{0}=70\,$km s-1 Mpc-1, $\Omega_{M}=0.3$,
and $\Omega_{\Lambda}=0.7$.

Figure 1: The recent SFHs (past 200 Myr) of three example simulated galaxies,
showing a range of $\eta$ values. The dark grey band encloses the previous
$10\,\mathrm{Myr}$, while the lighter band encloses the past
$100\;\mathrm{Myr}$; $\eta$ is the ratio of the SFRs averaged over those two
time periods. The black line shows an example of constant star formation with
$\eta\approx 0$, while the red line shows a ‘temporarily quenched’ case
($\eta=-0.64$). The blue line shows an example simulated galaxy observed
during a burst, where the star formation activity averaged over the previous
$10\;\mathrm{Myr}$ has increased compared to the average over the previous
$100\,\mathrm{Myr}$ (giving $\eta=0.88$).

##  2 Method

To provide the synthetic SEDs used for this test, we used four cosmological
zoom-in simulations of massive high-redshift galaxies from the FIRE
project.333http://fire.northwestern.edu The halos were first presented in
Feldmann et al. (2016). Subsequently, the halos were resimulated by Anglés-
Alcázar et al. (2017) using the FIRE-2 code (Hopkins et al., 2018); we use the
simulations from that work. The simulations use the code gizmo (Hopkins,
2015),444http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html with
hydrodynamics solved using the mesh-free Lagrangian Godunov ‘MFM’ method. The
simulations include cooling and heating from a meta-galactic background and
local stellar sources from $T\sim 10-10^{10}\,$K; star formation in locally
self-gravitating, dense, self-shielding molecular, Jeans-unstable gas; and
stellar feedback from OB & AGB mass-loss, SNe Ia & II, and multi-wavelength
photo-heating and radiation pressure, with inputs taken directly from stellar
evolution models. The FIRE physics, source code, and all numerical parameters
are exactly identical to those in Hopkins et al. (2018).

Snapshots of the halos’ central galaxies were taken from $1<z<8$ at intervals
of $15-25\,\mathrm{Myr}$. For each snapshot, we compute the true
$10\,\mathrm{Myr}$ and $100$ Myr-averaged SFRs,
$\overline{\psi}_{\mathrm{10\,Myr}}$ and
$\overline{\psi}_{\mathrm{100\,Myr}}$, respectively. Following Broussard et
al. (2019), we quantify the recent star forming activity using the ‘burst
indicator’
$\eta=\log_{10}(\mathrm{\overline{\psi}_{10\,\mathrm{Myr}}}/\mathrm{\overline{\psi}_{100\,\mathrm{Myr}}})$,
such that

$\eta\begin{cases}<0,&\text{galaxy has reduced SFR within the last
$10\,\mathrm{Myr}$};\\\ \approx 0,&\text{galaxy has had a constant SFR over
the last $100\,\mathrm{Myr}$};\\\ >0,&\text{galaxy has increased SFR in the
last $10\,\mathrm{Myr}$}.\\\ \end{cases}$

Figure 1 shows three example SFHs of snapshots used in this study with
$\eta\approx 0$ (black line), $\eta<0$ (red line) and $\eta>0$ (blue line).

Cochrane et al. (2019) calculated synthetic SEDs for these simulated galaxies
using the radiative transfer code SKIRT.555http://www.skirt.ugent.be/ Each
snapshot was ‘observed’ at 7 viewing angles spaced at 30∘ intervals, ranging
from aligned with the angular momentum vector to anti-aligned. This resulted
in 6,706 forward-modeled SEDs. See Cochrane et al. (2019); Cochrane et al.
(2022); Cochrane et al. (2023b); Cochrane et al. (2023a) and Parsotan et al.
(2021) for further details about the SKIRT calculations and other applications
of these and similar simulations/SEDs. We convolved these SEDs with the filter
response curves for 18 photometric bands spanning observed-frame wavelengths
$0.39<\lambda/\mu m<500$ (we refer the interested reader to H23 for further
details) assuming a signal-to-noise ratio of 5 in every band (following Smith
& Hayward, 2018). We fit the synthetic SEDs using two SED fitting codes,
Magphys (da Cunha et al., 2008, hereafter DC08) and Prospector (Leja et al.
2017, 2019, Johnson et al. 2021), fixing the redshift to the true value (i.e.
photometric redshift errors do not contribute to the bias demonstrated here).
We use the high-redshift version of Magphys (da Cunha et al., 2015) which
builds model UV–FIR SEDs by linking a stellar library containing 50,000 pre-
computed SFHs with those which satisfy the energy balance criterion among
another pre-computed library of 25,000 Charlot & Fall (2000) dust models; see
DC08 and H23 for detailed discussions. By calculating the $\chi^{2}$ goodness-
of-fit parameter between the observed photometry and every combination which
satisfies the energy balance, Magphys is able to marginalise over the prior
distribution to estimate posterior probability distributions for each property
in the model. Magphys assumes an initial mass function (IMF) from Chabrier
(2003), and stellar models from Bruzual & Charlot (2003) for metallicity
between 0.02 and 2 times solar.

For the Prospector fits, we followed the process described in Das et al., in
preparation. Prospector uses the Flexible Stellar Populations Synthesis (FSPS;
Conroy et al., 2009) code assuming a Kroupa (2001) IMF and solar metallicity.
Dust attenuation is modelled using the two-component Charlot & Fall (2000)
model and the Kriek & Conroy (2013) attenuation curve. The three-parameter
Draine & Li (2007, hereafter DL07) dust emission templates are used to model
the shape of the IR SED. The dynamic nested sampling library Dynesty (Speagle,
2020) is used to estimate Bayesian posteriors and evidences.

For both fitters, we consider an SED fit to be acceptable if the best fit
$\chi^{2}$ is below the 99 per cent confidence limit taken from standard
$\chi^{2}$ tables with the degrees of freedom calculated following Smith et
al. (2012, who calculated this only for Magphys; here we assume that it
applies equally well to the Prospector results, though this is perhaps
unlikely to be true given the differences between the two sets of models).
Throughout the subsequent analysis, we consider only those SEDs with
satisfactory fits.

Based on these analyses, we obtained four sets of estimated galaxy properties
for each snapshot. The Magphys values assume a fixed parameterisation of the
SFH based on an underlying continuous, delayed exponential model overlaid with
random starbursts (as detailed in DC08). In addition, we have three sets of
results from Prospector, assuming different SFH priors: “Dirichlet”,
“continuity”, and “bursty continuity”. The Dirichlet prior assumes that the
fractional specific star formation rate (sSFR) in each time bin follows a
Dirichlet distribution, with a concentration parameter, $\alpha$, that
controls the shape of the SFH (Leja et al., 2017, 2019), and which we set to
unity. The “continuity” prior directly fits for $\Delta\log(\mathrm{SFR})$
between neighbouring time bins using a Student’s $t$-distribution with scale
factor $\sigma=0.3$ and two degrees of freedom. This prior discourages abrupt
changes in SFR between adjacent bins but remains flexible enough to fit both
star-forming and quenched galaxies (Leja et al., 2019; Johnson et al., 2021).
Following Tacchella et al. (2022), we also consider the bursty continuity
prior, which is a modified form of the continuity prior with $\sigma=1$ and
$\nu=2$, shown to allow even greater degrees of variability in the shape of
the SFH. For each of these priors, we adopt the following binning scheme: the
first two and the last bins are kept the same for all sources, covering
$0<t_{l}<10$ Myrs, $10<t_{l}<100$ Myrs, and $0.85t_{z}<t_{l}<t_{z}$,
respectively. Here, $t_{z}$ represents the age of the Universe at the object’s
redshift, and $t_{l}$ is the lookback time. The remaining time between 100
Myrs and $0.85t_{z}$ is evenly spaced in logarithmic time. The number of bins
is chosen to ensure each bin has a width greater than
$0.02~{}\log_{10}\mathrm{Gyrs}$. For the lowest redshift snapshots in our
simulated data set, this corresponds to eight bins, while for the highest
redshifts, it is five.

For this study we are considering the fidelity of the $100\,\mathrm{Myr}$
averaged SFR inferred from the SED fitting of both Magphys and Prospector.
Following H23, we quantify the fidelity of the recovery of SFR, using the
residual in dex between the inferred value and the true value (from the
simulation) for the SFR, specifically:
$\Delta\log_{10}\psi=\log_{10}(\psi_{\mathrm{inferred}})-\log_{10}(\psi_{\mathrm{true}})$.

##  3 Results

Table 1 shows the fit success rate and the average SFR residual for each of
the four runs. It is immediately apparent that Prospector is able to produce
acceptable fits for a larger fraction of the snapshots ($>90$ per cent fit
success rate irrespective of the prior assumed) than Magphys (83 per cent).
This may reflect the effect that the Magphys pre-computed libraries – coupled
with the requirement that it consider only SFHs shorter than the age of the
Universe at a given redshift – result in the priors becoming increasingly
poorly-sampled at higher redshifts, and acceptable Magphys fits increasingly
hard to obtain at earlier epochs. Secondly, we note that the typical SFR
uncertainties for the Prospector runs shown in Table 1 are significantly
larger than those yielded by Magphys, and we will return to this point below.
These details notwithstanding, it is clear from Table 1 that there is no
significant difference in the overall $\Delta\log_{10}\psi$ once the
uncertainties are considered, meaning that both codes do a similarly
acceptable job of recovering the true SFR on average.

Table 1: SED fit success rate and typical SFR estimate fidelity, $\Delta\log_{10}(\psi)$, for the Magphys and Prospector fits to our synthetic observations. For the three sets of Prospector fits, the different SFH priors are indicated. Fitter/prior | % Fit success | $\Delta\log_{10}(\psi)$
---|---|---
Magphys | 83 | -0.11 $\pm$ 0.06
Prospector\- Bursty-continuity | 91 | -0.09 $\pm$ 0.24
Prospector\- Dirichlet ($\alpha=1$) | 92 | -0.30 $\pm$ 0.19
Prospector\- Continuity | 92 | -0.12 $\pm$ 0.12

We now investigate the fidelity of the Magphys and Prospector SFR estimates as
a function of the recent SFH. In Figure 2, we plot $\Delta\log(\psi)$ against
$\eta$, with the Magphys results shown in the left-hand panel, and the
Prospector results for the bursty continuity prior plotted in the right-hand
panel. The axes are chosen such that ‘starbursting’ snapshots appear to the
right of each plot, and ‘temp-quenched’ ones are on the left. The values for
individual SEDs are shown as blue points, with the red lines showing a running
average of $\Delta\log_{10}(\psi)$ for 100 data points. The red areas enclose
the 16th and 84th percentiles of the inferred SFR PDF using the same running
average, while the grey shading encloses the 16th and 84th percentiles of the
scatter within each corresponding subsample.

Figure 2: The LH panel shows the difference between the inferred and true SFRs
averaged over the last 100 Myr,
$\Delta\log_{10}\psi=\log_{10}(\psi_{\mathrm{inferred}})-\log_{10}(\psi_{\mathrm{true}})$,
as a function of the burst indicator
$\eta\equiv\log_{10}(\mathrm{\overline{\psi}_{10\,\mathrm{Myr}}}/\mathrm{\overline{\psi}_{100\,\mathrm{Myr}}})$.
The RH panel show the same obtained using Prospector with the Bursty
Continuity SFH prior (though similar results are obtained using different non-
parametric priors). In both panels the red line shows the mean residual
averaged over 100 data points with the red shaded area enclosing the averaged
16th and 84th percentile residuals from the PDF. The grey shaded area
indicates the 16th and 84th percentiles of the scatter in
$\Delta\log_{10}\psi$ within each averaged point, and the blue markers show
values for individual SEDs. The Magphys SFRs exhibit a systematic bias: the
SFR tends to be underestimated for ‘temp-quenched’ galaxies and overestimated
for galaxies that have experienced a burst within the past 10 Myr. This bias
is not apparent in the Prospector results.

Two effects are immediately apparent. First of all, there is a systematic
trend in the Magphys SFRs, such that while they are typically overestimated by
$\sim 0.4$ dex at $\eta>0.5$ (the ‘starbursting’ snapshots), the converse is
true at $\eta<-0.5$ (the ‘temp-quenched’ snapshots) where Magphys typically
underestimates the true SFRs by around 0.3 dex. This is clearly of concern,
since we know that both ‘starbursting’ and ‘temp-quenched’ galaxies exist in
the real Universe (e.g. Guo et al., 2016; Faisst et al., 2019; Broussard et
al., 2022; Looser et al., 2023) and this effect could – left unchecked –
result in a systematic error in the SFRs of galaxies that could be 1 dex in
magnitude in the worst cases. Secondly, and underscoring the severity of the
first issue, there is a clear tendency for the typical Magphys uncertainties
(shown as the red shading in the left panel) to underestimate the degree of
scatter visible in the data points (grey shading). This effect becomes
increasingly noticeable as $|\eta|$ increases.

The equivalent Prospector plot shows the SFR recovery –
$\Delta\log_{10}(\psi)$ – is approximately independent of $\eta$; the same
result is observed for the other Prospector priors with reasonable
uncertainties at $\eta_{\mathrm{true}}<0.5$ (where $\sim 90$ per cent of the
snapshots lie). At these $\eta$, the $\sim 0.3$ dex offset in the Prospector
SFRs relative to the ground truth is similar to that returned by Magphys at
$\eta<-0.5$; however, in Prospector the deviation from zero is not
statistically significant due to the more realistic error bars.

##  4 Discussion

Fundamental to SED fitting codes’ abilities to recover the true properties of
galaxies are (i) the wavelength sampling and (ii) the priors. The former issue
is closely related to the extent of our ability to estimate $\eta$ using SED
fitting (since different wavelength regimes have different ‘response
functions’ in terms of how the luminosity depends on the SFH; e.g. Kennicutt &
Evans, 2012; Sparre et al., 2017); we shall return to this topic in a future
work (Das et al., in preparation). In this work we focus solely on our ability
to recover 100 Myr-averaged SFRs using an example 18 bands of photometry
similar to that used in Smith et al. (2021) for the LOFAR deep fields. This
good wavelength sampling, in addition to the fact that we have not added noise
to photometry, means that this test represents a best-case scenario: the bias
that we have identified is fundamental, not due to noise or poor SED sampling,
and thus cannot be addressed with ‘better data’.

Figure 3: The distributions of $\eta$ values for 50,000 SFHs in the Magphys
high-$z$ libraries (in blue), overlaid with the ground truth values from the
simulations (hatched) as indicated in the legend. Note the significantly
larger $x$ range in this plot relative to Figure 2, and the logarithmic axes.

In the Magphys model, the priors are not under the control of the end user
since it uses pre-computed SFH libraries. The high-$z$ Magphys SFHs assume a
delayed exponential parametric form, with random bursts superposed in such a
manner that during the previous 2 Gyr, 75 per cent of SFHs include a burst
lasting between 30-300 Myr and forming between 0.1 and 100 times the stellar
mass formed by the underlying distribution (da Cunha et al., 2015). The $\eta$
distribution for the Magphys high-$z$ SFHs is shown by the blue histogram in
Figure 3, overlaid with the distribution of values for the true SFHs of the
simulated galaxies (hatched). It is clear that the libraries contain a
significantly broader range of $\eta$ values than the simulations, so the
effective prior on $\eta$ is not the source of the bias.

In Section 3 we also noted a difference between the typical uncertainties on
the Magphys SFR estimates and the observed scatter, indicating that the
uncertainties are increasingly underestimated as $|\eta|$ moves away from 0.
Figure 3 shows that that SFHs with such $\eta$ are far rarer in the Magphys
prior than those at $\eta\approx 0$, which may contribute to the
underestimated error bars (which are up to a factor $>5\times$ too small) in
the SFRs of these sources.

Prospector offers a solution to this problem, likely as a result of the more
flexible non-parametric SFH employed, and the improved ability to explore the
parameter space using dynamic nested sampling. Even so, our results
demonstrate that the choice of prior can impact the error estimates associated
with Prospector SFRs for the most bursty galaxies (those with highest $\eta$),
which is unsurprising given e.g. that the continuity prior tends to down-
weight rapid changes in the SFHs.

##  5 Conclusions

We have used four cosmological zoom-in simulations spanning $1<z<8$, together
with the radiative transfer code SKIRT, to generate 6,706 synthetic galaxy
SEDs. We attempted to recover the true properties of synthetic observations
based on these simulations using Magphys and Prospector. Although there is
little difference between the fidelity of SFR estimates inferred using the two
codes when averaging over the whole sample, we find that the accuracy of the
Magphys-inferred SFR is strongly dependent on the recent SFH. This trend is
sufficiently large that the Magphys 100 Myr-averaged SFRs of individual
galaxies could be overestimated (underestimated) by as much as 1 dex in the
worst cases for galaxies experiencing starbursts (temporarily quenched
periods). The trend is not evident in the 100 Myr SFRs from Prospector, which
we speculate is due to the non-parametric SFHs assumed and possibly the
dynamic nested sampling method used to better explore the parameter space.
Given this bias, we caution against employing parametric SFHs when fitting the
SEDs of observed galaxies that are likely to have bursty SFHs, such as high-
redshift galaxies observed with JWST.

## Acknowledgements

PH and DJBS would like to thank the Flatiron Institute for their hospitality.
DJBS acknowledges support from the UK Science and Technology Facilities
Council (STFC) under grant ST/V000624/1. This research has made use of NASA’s
Astrophysics Data System Bibliographic Services.

## Data Availability

The data presented in this paper will be shared on reasonable request to the
first author.

## References

  * Anglés-Alcázar et al. (2017) Anglés-Alcázar D., Faucher-Giguère C.-A., Quataert E., Hopkins P. F., Feldmann R., Torrey P., Wetzel A., Kereš D., 2017, MNRAS, 472, L109
  * Baes (2019) Baes M., 2019, Proceedings of the International Astronomical Union, 15, 26
  * Baldwin et al. (2018) Baldwin C., McDermid R. M., Kuntschner H., Maraston C., Conroy C., 2018, MNRAS, 473, 4698
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57
  * Best et al. (2023) Best P. N., et al., 2023, MNRAS, 523, 1729
  * Broussard et al. (2019) Broussard A., et al., 2019, ApJ, 873, 74
  * Broussard et al. (2022) Broussard A., Gawiser E., Iyer K., 2022, ApJ, 939, 35
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Calzetti (1997) Calzetti D., 1997, AJ, 113, 162
  * Carnall et al. (2019) Carnall A. C., Leja J., Johnson B. D., McLure R. J., Dunlop J. S., Conroy C., 2019, ApJ, 873, 44
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Charlot & Fall (2000) Charlot S., Fall S. M., 2000, ApJ, 539, 718
  * Cochrane et al. (2019) Cochrane R. K., et al., 2019, MNRAS, 488, 1779
  * Cochrane et al. (2021) Cochrane R. K., et al., 2021, MNRAS, 503, 2622
  * Cochrane et al. (2022) Cochrane R. K., Hayward C. C., Anglés-Alcázar D., 2022, ApJ, 939, L27
  * Cochrane et al. (2023a) Cochrane R. K., Anglés-Alcázar D., Cullen F., Hayward C. C., 2023a, arXiv e-prints, p. arXiv:2310.08829
  * Cochrane et al. (2023b) Cochrane R. K., Hayward C. C., Anglés-Alcázar D., Somerville R. S., 2023b, MNRAS, 518, 5522
  * Conroy (2013) Conroy C., 2013, ARA&A, 51, 393
  * Conroy et al. (2009) Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486
  * da Cunha et al. (2008) da Cunha E., Charlot S., Elbaz D., 2008, MNRAS, 388, 1595
  * da Cunha et al. (2015) da Cunha E., et al., 2015, ApJ, 806, 110
  * Draine & Li (2007) Draine B. T., Li A., 2007, ApJ, 657, 810
  * Dudzevičiūtė et al. (2020) Dudzevičiūtė U., et al., 2020, MNRAS, 494, 3828
  * Dye (2008) Dye S., 2008, MNRAS, 389, 1293
  * Faisst et al. (2019) Faisst A. L., Capak P. L., Emami N., Tacchella S., Larson K. L., 2019, ApJ, 884, 133
  * Feldmann et al. (2016) Feldmann R., Hopkins P. F., Quataert E., Faucher-Giguère C.-A., Kereš D., 2016, MNRAS, 458, L14
  * Flores Velázquez et al. (2021) Flores Velázquez J. A., et al., 2021, MNRAS, 501, 4812
  * Guo et al. (2016) Guo Y., et al., 2016, ApJ, 833, 37
  * Haskell et al. (2023) Haskell P., Smith D. J. B., Cochrane R. K., Hayward C. C., Anglés-Alcázar D., 2023, MNRAS, 525, 1535
  * Hayward & Hopkins (2017) Hayward C. C., Hopkins P. F., 2017, MNRAS, 465, 1682
  * Hayward & Smith (2015) Hayward C. C., Smith D. J. B., 2015, MNRAS, 446, 1512
  * Hodge et al. (2016) Hodge J. A., et al., 2016, ApJ, 833, 103
  * Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  * Hopkins et al. (2023) Hopkins P. F., et al., 2023, MNRAS, 525, 2241
  * Hunt et al. (2019) Hunt L. K., et al., 2019, A&A, 621, A51
  * Iyer et al. (2020) Iyer K. G., et al., 2020, MNRAS, 498, 430
  * Johnson et al. (2021) Johnson B. D., Leja J., Conroy C., Speagle J. S., 2021, ApJS, 254, 22
  * Kennicutt & Evans (2012) Kennicutt R. C., Evans N. J., 2012, ARA&A, 50, 531
  * Kokorev et al. (2021) Kokorev V. I., et al., 2021, ApJ, 921, 40
  * Kriek & Conroy (2013) Kriek M., Conroy C., 2013, ApJ, 775, L16
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kroupa & Jerabkova (2019) Kroupa P., Jerabkova T., 2019, Nature Astronomy, 3, 482
  * Lee (2010) Lee S.-K. J., 2010, PhD thesis, Johns Hopkins University, Maryland
  * Lee et al. (2009) Lee S.-K., Idzi R., Ferguson H. C., Somerville R. S., Wiklind T., Giavalisco M., 2009, ApJS, 184, 100
  * Leja et al. (2017) Leja J., Johnson B. D., Conroy C., van Dokkum P. G., Byler N., 2017, ApJ, 837, 170
  * Leja et al. (2019) Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2019, ApJ, 876, 3
  * Liu (2020) Liu D., 2020, michi2: SED and SLED fitting tool, Astrophysics Source Code Library, record ascl:2005.002 (ascl:2005.002)
  * Looser et al. (2023) Looser T. J., et al., 2023, arXiv e-prints, p. arXiv:2306.02470
  * Lower et al. (2020) Lower S., Narayanan D., Leja J., Johnson B. D., Conroy C., Davé R., 2020, ApJ, 904, 33
  * Madau et al. (1996) Madau P., Ferguson H. C., Dickinson M. E., Giavalisco M., Steidel C. C., Fruchter A., 1996, MNRAS, 283, 1388
  * Maraston (2005) Maraston C., 2005, MNRAS, 362, 799
  * Mobasher et al. (2015) Mobasher B., et al., 2015, ApJ, 808, 101
  * Noeske et al. (2007) Noeske K. G., et al., 2007, ApJ, 660, L43
  * Noll et al. (2009) Noll S., Burgarella D., Giovannoli E., Buat V., Marcillac D., Muñoz-Mateos J. C., 2009, A&A, 507, 1793
  * Pacifici et al. (2023) Pacifici C., et al., 2023, ApJ, 944, 141
  * Parsotan et al. (2021) Parsotan T., Cochrane R. K., Hayward C. C., Anglés-Alcázar D., Feldmann R., Faucher-Giguère C. A., Wellons S., Hopkins P. F., 2021, MNRAS, 501, 1591
  * Pforr et al. (2012) Pforr J., Maraston C., Tonini C., 2012, MNRAS, 422, 3285
  * Rujopakarn et al. (2019) Rujopakarn W., et al., 2019, ApJ, 882, 107
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Simha et al. (2014) Simha V., Weinberg D. H., Conroy C., Dave R., Fardal M., Katz N., Oppenheimer B. D., 2014, arXiv e-prints, p. arXiv:1404.0402
  * Smith & Hayward (2015) Smith D. J. B., Hayward C. C., 2015, MNRAS, 453, 1597
  * Smith & Hayward (2018) Smith D. J., Hayward C. C., 2018, MNRAS, 476, 1705
  * Smith et al. (2012) Smith D. J. B., et al., 2012, MNRAS, 427, 703
  * Smith et al. (2021) Smith D. J. B., et al., 2021, A&A, 648, A6
  * Sparre et al. (2017) Sparre M., Hayward C. C., Feldmann R., Faucher-Giguère C.-A., Muratov A. L., Kereš D., Hopkins P. F., 2017, MNRAS, 466, 88
  * Speagle (2020) Speagle J. S., 2020, MNRAS, 493, 3132
  * Tacchella et al. (2022) Tacchella S., et al., 2022, ApJ, 927, 170
  * Walcher et al. (2011) Walcher J., Groves B., Budavári T., Dale D., 2011, Ap&SS, 331, 1
  * Wang et al. (2022) Wang Y., et al., 2022, MNRAS, 515, 3249
  * Wuyts et al. (2009) Wuyts S., Franx M., Cox T. J., Hernquist L., Hopkins P. F., Robertson B. E., van Dokkum P. G., 2009, ApJ, 696, 348