Further author information: (Send correspondence to J.F.)
J.F.: E-mail<EMAIL_ADDRESS>
# Closed-Loop Until Further Notice: Comparing Predictive Control Methods in
Closed-Loop

J. Fowler  Department of Astronomy & Astrophysics, University of California,
Santa Cruz, CA, USA  M.A.M. van Kooten  National Research Council Canada,
Herzberg Astronomy and Astrophysics Research Center, Victoria, Canada  R.
Jensen-Clem  Department of Astronomy & Astrophysics, University of California,
Santa Cruz, CA, USA

###### Abstract

For future extremely large telescopes, error in extreme adaptive optics
systems at small angular separations will be highly impacted by the lag time
of the correction, which is typically on millisecond timescales; one solution
is to apply a predictive correction to catch up with the system delay.
Predictive control leads to significant RMS error reductions in simulation (on
the order of 5-10x improvement in RMS error compared with a standard integral
controller), but shows only modest improvement on-sky (less than 2x in RMS
error). This performance limitation is likely impacted by elements of pseudo
open loop (POL) reconstruction, which requires assumptions about the response
of the deformable mirror and accuracy of the wavefront measurements that are
difficult to verify in practice. In this work, we explore a closed-loop method
for data-driven prediction using a reformulated empirical orthogonal functions
(EOF). We examine the performance of the open and closed-loop methods in
simulation on perfect systems and systems with an inaccurate understanding of
the DM response.

###### keywords:

predictive control, extreme adaptive optics (exAO), empirical orthogonal
functions (EOF)

##  1 INTRODUCTION

Dessenne, 1998 [3] was one of the first papers outlining predictive control
for adaptive optics (AO) systems. This predictive method was demonstrated on-
sky in 1999 [2] at the 1.52 meter observatory in Haute Provence; they cite a
relative Strehl increase over a classic integrator of 30$\%$ (a max Strehl
performance of $\sim 14\%$) at their visible central wavelength of 650 nm.
Since then, a plethora of predictive control methods have come into the
literature, including linear estimators [5, 6, 9, 22], linear quadratic
Gaussian controllers [13, 19, 18, 14], model-based updates to a linear
quadratic Gaussian [16, 17, 4], subspace control methods [8], and non-linear
neural network solvers [11, 23, 10, 20, 12, 7]. Despite nearly 25 years of
predictive methods, Empirical Orthogonal Functions (EOF)[5] is the only method
to have been run on-sky as the controller for all spatial frequencies on an
8-10 meter class telescope; demonstrated on Subaru/SCExAO[6] and on Keck/NIRC2
[22].

As opposed to more classic methods (e.g., an integral controller or a linear
quadratic Gaussian controller) a classic EOF [5] learns a linear relationship
in the evolution of the full scale of input turbulence, which means that EOF
trains on (and predicts) the open loop wavefront. Because no single conjugate
AO system runs in open loop, in practice these methods train on pseudo-open
loop (POL) data, where deformable mirror (DM) commands are added to wavefront
sensor (WFS) measurements to reconstruct the full state of turbulence.
However, this leaves room for error in the reconstruction. Any spatially or
temporally evolving mismatches in the calibration (e.g., non-linearities in
the sensing and actuator pokes or DM to WFS misregistrations) that are not as
impactful in an integral closed-loop controller could lead to incorrect
performance prediction in open loop.

However, upon closer inspection of Dessenne’s adaptive predictor [3] and
Guyon’s EOF [5], it becomes clear that they are both methods that build
predictive controllers from the same pieces of information. The major
departure between the two methods is that Guyon’s EOF runs in open loop, and
reconstructs full turbulence states, and Dessenne’s preditive controller runs
in closed-loop, and encapsulates the evolution of both the DM commands and WFS
measurements. With this work, we aim to revisit Dessenne’s method, update it
for comparison with Guyon’s EOF, perform preliminary simulations of its
feasibility within the context of a modern system, and ultimately begin to
compare what is essentially an open and closed-loop implementation of
empirical orthogonal functions. From hence forth, we will refer to both
methods as EOF, and make the distinction between open and closed-loop
implementations.

In Section 2.1, we describe the pseudo-open loop EOF method, in Section 2.2 we
describe the closed-loop update for this work, and in Section 2.3 we discuss
training conditions for the closed-loop update. In Section 3 we present the
results of preliminary performance simulations and the impact of model-
mismatch errors on that performance.

##  2 EMPIRICAL ORTHOGONAL FUNCTIONS

Empirical orthogonal functions predicts a future state of a wavefront by
building a predictive filter that linearly combines previous states of the
wavefront. Figure 1 shows a visual representation. In the following sections
we describe the math that represents this process for open and closed-loop
implementations.

---
Figure 1:  The future state of the wavefront is predicted as a linear
combination of previous states, $w$ at some time $n$. The collection of
previous states is called the history vector, $\mathbf{h}$. Phase screens
evolve in time (in sea foam) and a weight is applied to each one (the
predictive filter $\mathbf{F}=[f_{1},f_{2},...]$, where each $f_{n}$ contains
m weights for m modes) to estimate the final wavefront prediction in teal.
Note that mathematically a history vector is a single flat vector that
contains the appended information of all 3 screens to predict a final vector
that contains a flattened version of the predicted wavefront.

###  2.1 Previous Open Loop Implementations

Starting with Guyon, 2017 [5] and a follow up from Jensen-Clem in 2019 [9], we
outline our implementation of EOF in (pseudo-)open loop. Given some state of
the wavefront (i.e., full scale of the uncorrected turbulence as considered in
pseudo open loop), $w$, with m variables representing wavefront sensor
measurements (which for a zonal approach corresponds to the number of
deformable mirror actuators) and n associated frames, we build a history
vector for each subaperture:

$\mathbf{h}_{m}(t)=\begin{bmatrix}w(t)\\\ w(t-dt)\\\
w(t-(n-1)dt)\end{bmatrix}$ (1)

We build a predictive filter for each mode Fm that will predict the full phase
at a given point with:

$\mathbf{F}_{m}\mathbf{h}_{m}(t)=w(t+dt)$ (2)

Both the filter and history vector can be written to include information
across multiple modes, i.e.,
$\mathbf{h}=[w_{0}(t),w_{1}(t)...w_{0}(t-dt),w_{1}(t-dt)...]$, but for ease of
comparison with the closed-loop implementation we leave these as single
wavefront sensor measurement/DM actuator filters and predictions.

To find the predictive filter $\mathbf{F}$ we minimize an error term that
consists of the difference between the output predicted wavefront (in DM
space) and the true phase at that time. We collect training data $\mathbf{D}$,
which contains history vectors and $\mathbf{P}$ which holds the true state (or
“future”) for each history vector. (I.e., mapping $t$ to its future state at
time $t+dt$.)

$\textrm{min}_{\mathbf{F}}||\mathbf{D}^{T}\mathbf{F}^{T}-\mathbf{P}^{T}||^{2}$
(3)

Solving equation (3) requires a pseudo-inverse; Guyon, 2017 [5] solved this
problem with an SVD (singular value decomposition) inversion, but we use a
least-squares inversion [9], with regularization constant $\alpha$ ($\alpha$
may be set to $1$ for simulations, but is found empirically on-sky [22].)

$\displaystyle\mathbf{F}$ $\displaystyle=$
$\displaystyle((\mathbf{D}^{T})^{\dagger}\mathbf{P}^{T})^{T}$ (4)
$\displaystyle\mathbf{F}$ $\displaystyle=$
$\displaystyle\mathbf{P}\mathbf{D}^{T}(\mathbf{D}\mathbf{D}^{T}+\alpha\mathbf{I})^{-1}$
(5)

Finally, our predictive filter $\mathbf{F}$ holds a coefficient for each
previous state, as expressed by a pseudo open loop (POL) reconstruction of the
telemetry projected into DM space.

###  2.2 Our Development of a Closed-Loop Update

Inspired by Dessenne, 1998 [3] and Haffert, 2021 [8], we explore a closed-loop
reformulation of the classic empirical orthogonal functions [5]. The fullest
realization of this method could:

  1. 1.

Improve controller stability by running in closed-loop.

  2. 2.

Avoid error introduced by non-linear wavefront sensing or DM model mismatch
when wavefront sensor residuals are converted to psudeo open loop.

  3. 3.

Track drifts that may impact the DM and wavefront sensor separately by
allowing each to evolve with a different set of coefficients.

  4. 4.

Do robust timekeeping with information into the system.

Following [3], an updated history vector (we refer to this as
$\boldsymbol{\phi}$ to distinguish from the open loop history vector, as in
standard in control theory conventions), contains both wavefront sensor
residuals $\epsilon(t)$ and DM commands $y(t)$. (While the original work opted
to use KL-modes, we work in the DM zonal basis, i.e., one point of information
per DM actuator/WFS subaperture.) We consider only a single mode at a time:

$\boldsymbol{\phi}(n)=\begin{bmatrix}y(n-1)\\\ y(n-2)\\\ ...\\\ y(n-p)\\\
\epsilon(n-2)\\\ \epsilon(n-3)\\\ ...\\\ \epsilon(n-p-1)\end{bmatrix}$ (6)

In this way, the new history vector and the corresponding predictive filter
have twice as many values per mode pair as previous open loop derivations.
However, the output of the filter applied to the history vector
$\boldsymbol{\theta}\cdot\boldsymbol{\phi}$ predicts the same piece of
information as an open loop implementation: the full turbulence in DM space,
which is the DM command needed at a given iteration, in our notation $y(n)$.

Dessenne’s original formulation suggests reconstructing open loop turbulence
from the full transfer function of the AO control loop to train the predictive
filter, and solving for a steady state recursive least squares solution [3].
In this simulation, we instead apply the same minimization technique used in
the open loop implementation. Future work, inspired by that of van Kooten,
2019 [21] will explore a recursive least square implementation that updates
the filter with each control iteration. We build a new minimization problem:

$\textrm{min}_{\boldsymbol{\theta}}||\mathbf{D}^{T}\boldsymbol{\theta}^{T}-\mathbf{P}^{T}||^{2}$
(7)

where $\mathbf{D}$ contains collections of the history vector
$\boldsymbol{\phi}$, and $\mathbf{P}$ contains the future state of the full
turbulence.

In practice, the distinction between the open and closed-loop implementation
is twofold: (1) the ability of the correction to apply different coefficients
for the DM commands and the wavefront sensor information (whereas POL is one
command applied to a summed value) and (2) the ability to robustly encapsulate
a time delay into the reconstruction. POL reconstruction adds wavefront sensor
measurements to DM commands with a single static delay, while closed-loop
methods can account for when each piece of information enters the system (with
two-step delays or delays with non-uniform steps). See Appendix A for
additional information on how timekeeping impacts the construction of
$\boldsymbol{\phi}$ and $\mathbf{P}$.

Assuming a classic two-step delay and $p$ associated frames of information in
a history vector, the weighted estimate for a single coefficient for the
closed-loop formulation ($y_{CL}$) vs. the original open loop formulation
($y_{OL}$) takes the form:

$\displaystyle y_{CL}(n)$ $\displaystyle=$ $\displaystyle
b_{1}y(n-1)+...+b_{p}y(n-p)+a_{0}\epsilon(n-2)+...+a_{p-1}\epsilon(n-p-1)$ (8)
$\displaystyle y_{OL}(n)$ $\displaystyle=$ $\displaystyle
c_{0}[y(n-1)+\epsilon(n-1)]+...+c_{p}[y(n-p-1)+\epsilon(n-p-1)]$ (9)

###  2.3 Selecting a Truth Condition for Training Data

The truth condition P, against which the predictive filter trains, must be the
full uncorrected state of the wavefront (i.e., open loop turbulence) at a
given iteration. The original work from Dessenne [3] estimated the open loop
wavefront from the transfer function of the AO control loop, essentially a
higher fidelity pseudo open loop reconstruction than that of Guyon, 2017 [5].
For these simulations, we give the predictor perfect knowledge, using the full
scale of turbulence from simulation. We note that this training condition is
not realistic to on-sky operation, but acts as a laboratory to test the
perfect performance of a closed-loop implementation, and future work will
explore more realistic P generation. Figure 2 shows the example data used as
the truth condition for training the closed-loop implementation.

---
Figure 2:  Left: Piston, tip, and tilt subtracted simulated phase screens from
an 8 meter telescope with a single wind layer. Right: Table of parameters for
generating the turbulence and AO simulations.

##  3 PRELIMINARY IMPLEMENTATION OF CLOSED-LOOP PREDICTIVE CONTROL

We simulate atmospheric phase screens and an idealized AO control system with
HCIPy[15]. We estimate the root-mean-square (RMS) residual error across the
pupil for 1) the full state of the uncorrected turbulence, 2) a pseduo-
integrator, in which we apply a perfect correction with a 2 ms time delay (for
Keck, the time-delay is $\sim$1.5 ms[1]), 3) open loop EOF (well optimized in
previous work [9, 4]), and 4) a preliminary implementation of closed-loop EOF.
These simulations use perfect wavefront sensing and correction, and are meant
to provide comparative estimates only for the bandwidth error. The parameters
for the simulation are shown in Figure 2.

---
Figure 3:  Left: Two predictive methods alongside a psuedo-integrator (a
perfect correction 2 timesteps behind), all applied with no measurement or
fitting error (i.e., perfect wavefront sensing and control) running at 1 kHz.
Full uncorrected turbulence RMS median: 1534.21 nm, quasi-integrator: 24.30
nm, closed-loop predictor: 9.49 nm, open loop predictor: 1.03 nm. Right:
Impacts of a DM-model mismatch. Closed-loop EOF correction goes from a median
of 9.49 (for the nominal run) to 16.66 nm RMS when a DM model error is
introduced, but still shows better performance than the pseudo-integrator. The
open loop EOF error increases by an order of magnitude and is not displayed on
the figure for scale purposes.

###  3.1 Performance in an Idealized AO System

Figure 3 (left) shows the results of our initial implementation; we find a
factor of improvement in RMS error of $\sim$ 2.5 over a standard integrator
from closed-loop EOF, under-performing compared to the standard psuedo open
loop. In practice, we should be able to recreate similar (if not better)
performance between a closed and open loop implementation, however closed-loop
minimization problem contains twice as many variables (treating the WFS and DM
independently), twice as many regressors, and will require more training data
to converge to an optimal solution, as shown in van Kooten, 2019 [21]. Future
work will explore optimizing filter length and training length for a closed-
loop implementation for a more fair comparison with open loop.

###  3.2 Robustness to Model Errors

Figure 3 (right) shows the performance of both predictors on a second set of
data, in which we simulate a DM model mismatch in our control system – adding
a static factor of 2 between what the DM expects to apply and actually
applies. We see that the performance of closed-loop EOF is less sensitive to
errors in the system model. While the closed-loop performance drops when a DM
model issue is introduced, the closed-loop predictor still outperforms a
typical integrator. However, introducing the same issue into a open loop EOF
increases the residual error by more than an order of magnitude. We note that
in theory a closed-loop EOF should exactly learn and reconstruct the model-
error, and future work to optimize the filter and training length will likely
make a closed-loop EOF even more robust to model errors.

##  4 CONCLUSION

In conclusion, we revisit the work of Dessenne, 1998 [3] to compare open and
closed-loop implementations of empirical orthogonal function (EOF) in
simulation. A preliminary simulation of closed-loop EOF does not perform as
well as an optimized implementation of open loop EOF, but still provides
improvement over a classic integrator. A closed-loop implementation also
proves to be more robust to the introduction of DM model errors.

We speculate that future optimization of the closed-loop predictive filter
(e.g., exploring history vector length, training data length, etc) will likely
close the gap between the open and closed-loop implementations and provide a
controller that is even more robust to model errors. We also note that for our
preliminary simulations we used perfect knowledge of the system to train our
closed-loop predictive filter, which is not realistic for on-sky
implementation; if performance comparisons prove promising, a more realistic
training method could be devised, for example using the transfer function as
in Desenne’s original work [3] or training on open and closed-loop data at the
beginning of the night. Future work could also examine the comparative
stability of open and closed-loop methods, as well as the impact of a more
robust consideration of time-delay.

We revisit this closed-loop predictive controller not only as an exploratory
method for AO control, but also as a laboratory to explore model-mismatch and
improve performance of the open loop implementations that are operating on-
sky. For extreme adaptive optics on extremely large telescopes, novel ways to
account for bandwidth error are worth pursuing.

##  Appendix A TIME REFERENCE FRAMES AND TURBULENCE RECONSTRUCTION

One intriguing issue when comparing and reworking control methods is a robust
understanding of time delays over the course of a control loop and how that
plays out with different methods. Pseudo-open loop implementations gloss over
this issue, by recreating a single point of information at a single point in
time, but closed-loop implementations provide the opportunity to represent
each piece of information at the time it enters the system, some work even
accounts for fractional time-delays [17]. The following is a minimal proof of
information flow in a control loop, switching between time reference frames.

---
Figure 4:  Control diagram of a classic two-step delay. The wavefront ($s$)
comes into the system, and first meets a sum junction, where we apply
deformable mirror (DM) commands ($y$), making it a closed-loop system. At this
point, our residual is $\epsilon=s+y$. That signal is sensed by a wavefront
sensor; that sensed signal is $\epsilon^{\prime}=\textrm{WF}(\epsilon)$, where
WF represents a functional form of how the wavefront sensor interprets the
signal. Finally, this is fed into a computer that will calculate and apply a
correction $y=DM(\epsilon^{\prime})$, where DM represents a functional form of
how that correction is applied and calculated based on $\epsilon^{\prime}$.

We consider this system in two reference frames: (1) the time represented by
the physical information flowing through the system (i.e., what time the
turbulence the system is analyzing occurred); this is physically intuitive for
timekeeping and residual comparison and we call it a physical clock and (2)
the time represented by the control system and when information could be
sampled from various sensors or correctors, which we call a control clock.

Dessenne 1998 [3], uses a control clock framework. If, for example, we wanted
to associated 3 states in time, we would build a history vector (used to
predict the state at iteration n) of the form:

$\boldsymbol{\phi}(n)=\begin{bmatrix}y(n-1)\\\ y(n-2)\\\ y(n-3)\\\
\epsilon^{\prime}(n-2)\\\ \epsilon^{\prime}(n-3)\\\
\epsilon^{\prime}(n-4)\end{bmatrix}$ (10)

The goal of this appendix is to show that these time steps from the
perspective of the control clock associate logical pieces of information from
the perspective of a physical clock. If we consider information moving through
a system with a two step delay, we would see the chain of control events
outlined in Table 1

Table 1: How information moves through the control loop, starting from the loop turning on. With a step of delay between each component, we at first have no information sensed and turned into a DM command.  $s$ | $\epsilon$ | $\epsilon^{\prime}$ | $y$
---|---|---|---
$s_{1}$ | $\epsilon_{1}=s_{1}$ | - | -
$s_{2}$ | $\epsilon_{2}=s_{2}$ | $\epsilon^{\prime}_{2}=\textrm{WF}(\epsilon_{1})$ | -
$s_{3}$ | $\epsilon_{3}=s_{3}+y_{3}$ | $\epsilon^{\prime}_{3}=\textrm{WF}(\epsilon_{2})$ | $y_{3}=\textrm{DM}(\epsilon^{\prime}_{2})$
$s_{4}$ | $\epsilon_{4}=s_{4}+y_{4}$ | $\epsilon^{\prime}_{4}=\textrm{WF}(\epsilon_{3})$ | $y_{4}=\textrm{DM}(\epsilon^{\prime}_{3})$
$s_{5}$ | $\epsilon_{5}=s_{5}+y_{5}$ | $\epsilon^{\prime}_{5}=\textrm{WF}(\epsilon_{4})$ | $y_{5}=\textrm{DM}(\epsilon^{\prime}_{4})$
$s_{6}$ | $\epsilon_{6}=s_{6}+y_{6}$ | $\epsilon^{\prime}_{6}=\textrm{WF}(\epsilon_{5})$ | $y_{6}=\textrm{DM}(\epsilon^{\prime}_{5})$
$s_{7}$ | $\epsilon_{7}=s_{7}+y_{7}$ | $\epsilon^{\prime}_{7}=\textrm{WF}(\epsilon_{6})$ | $y_{7}=\textrm{DM}(\epsilon^{\prime}_{6})$
$s_{8}$ | $\epsilon_{8}=s_{8}+y_{8}$ | $\epsilon^{\prime}_{8}=\textrm{WF}(\epsilon_{7})$ | $y_{8}=\textrm{DM}(\epsilon^{\prime}_{7})$

If we now consider the history vector $\boldsymbol{\phi}$, we can rewrite it
to prediction some iteration $n=7$

$\boldsymbol{\phi}(7)=\begin{bmatrix}y(n-1)=y_{6}=\textrm{DM}(\epsilon_{5}^{\prime})=\textrm{DM(WF}(\epsilon_{4}))\\\
y(n-2)=y_{5}=\textrm{DM}(\epsilon_{4}^{\prime})=\textrm{DM(WF}(\epsilon_{3}))\\\
y(n-3)=y_{4}=\textrm{DM}(\epsilon_{3}^{\prime})=\textrm{DM(WF}(\epsilon_{2}))\\\
\epsilon^{\prime}(n-2)=\epsilon_{5}^{\prime}=\textrm{WF}(\epsilon_{4})\\\
\epsilon^{\prime}(n-3)=\epsilon_{4}^{\prime}=\textrm{WF}(\epsilon_{3})\\\
\epsilon^{\prime}(n-4)=\epsilon_{3}^{\prime}=\textrm{WF}(\epsilon_{2})\end{bmatrix}$
(11)

Notice that though the control clock indexing appears to introduce a one step
offset, the most expanded part of the expression for each iteration shows that
the elements on the history vector depend on the the same physically timed
pieces of information.

It should be noted that a classic two-step delay is not actually physically
reminiscent of most systems. In Cetre, 2018[1] they found that for the Keck
pyramid WFS RTC, the time delay over the entire correction (including
wavefront sensing time, calculations, and DM latency after correction is
applied) takes $\sim$1.5 ms, where $\sim$1 ms is the wavefront sensing and
readout time, and $\sim$0.5 ms is the calculation and hardware latency time.
For this simulation work we have opted to only update the DM once per new
piece of wavefront sensor information, essentially forcing this to be a
classic two-step delay. However, future work will explore the efficacy of
updating the DM more frequently, as we could project forward the correction
easily and provide two DM updates per wavefront sensor readout.

## References

  * [1] Sylvain Cetre et al.  “A near-infrared pyramid wavefront sensor for Keck adaptive optics: real-time controller”  In _Adaptive Optics Systems VI_ 10703, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 2018, pp. 1070339  DOI: 10.1117/12.2311781
  * [2] C. Dessenne, P.-Y. Madec and G. Rousset  “Sky implementationof modal predictive control in adaptive optics”  In _Optics Letters_ 24.5, 1999, pp. 339–341  DOI: 10.1364/OL.24.000339
  * [3] Caroline Dessenne, Pierre-Yves Madec and Gérard Rousset  “Optimization of a Predictive Controller for Closed-Loop Adaptive Optics”  In _Applied Optics_ 37.21, 1998, pp. 4623–4633  DOI: 10.1364/AO.37.004623
  * [4] J. Fowler, M… Van Kooten and R. Jensen-Clem  “Battle of the predictive wavefront controls: comparing data and model-driven predictive control for high contrast imaging”  In _Adaptive Optics Systems VIII_ 12185, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 2022, pp. 1218582  DOI: 10.1117/12.2629521
  * [5] Olivier Guyon and Jared Males  “Adaptive Optics Predictive Control with Empirical Orthogonal Functions (EOFs)”  In _arXiv e-prints_ , 2017, pp. arXiv:1707.00570  DOI: 10.48550/arXiv.1707.00570
  * [6] Olivier Guyon et al.  “The compute and control for adaptive optics (CACAO) real-time control software package”  In _Adaptive Optics Systems VI_ 10703, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 2018, pp. 107031E  DOI: 10.1117/12.2314315
  * [7] Rehan Hafeez et al.  “Forecasting wavefront corrections in an adaptive optics system”  In _Journal of Astronomical Telescopes, Instruments, and Systems_ 8, 2022, pp. 029003  DOI: 10.1117/1.JATIS.8.2.029003
  * [8] Sebastiaan Y. Haffert et al.  “Data-driven subspace predictive control of adaptive optics for high-contrast imaging”  In _Journal of Astronomical Telescopes, Instruments, and Systems_ 7, 2021, pp. 029001  DOI: 10.1117/1.JATIS.7.2.029001
  * [9] Rebecca Jensen-Clem et al.  “Demonstrating predictive wavefront control with the Keck II near-infrared pyramid wavefront sensor”  In _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ 11117, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 2019, pp. 111170W  DOI: 10.1117/12.2529687
  * [10] Rico Landman, Sebastiaan Y. Haffert, Vikram M. Radhakrishnan and Christoph U. Keller  “Self-optimizing adaptive optics control with reinforcement learning for high-contrast imaging”  In _Journal of Astronomical Telescopes, Instruments, and Systems_ 7, 2021, pp. 039002  DOI: 10.1117/1.JATIS.7.3.039002
  * [11] Xuewen Liu et al.  “Wavefront prediction using artificial neural networks for open-loop adaptive optics”  In _Monthly Notices of the Royal Astronomical Society_ 496.1, 2020, pp. 456–464  DOI: 10.1093/mnras/staa1558
  * [12] J. Nousiainen et al.  “Toward on-sky adaptive optics control using reinforcement learning. Model-based policy optimization for adaptive optics”  In _Astronomy & Astrophysics_ 664, 2022, pp. A71  DOI: 10.1051/0004-6361/202243311
  * [13] Randall N. Paschall and David J. Anderson  “Linear quadratic Gaussian control of a deformable mirror adaptive optics system with time-delayed measurements”  In _Applied Optics_ 32.31, 1993, pp. 6347–6358  DOI: 10.1364/AO.32.006347
  * [14] C. Petit et al.  “Integration, tests and laboratory performance of SAXO, the VLT-SPHERE extreme AO system”  In _Second International Conference on Adaptive Optics for Extremely Large Telescopes. Online at ¡A href=“http://ao4elt2.lesia.obspm.fr”¿http://ao4elt2.lesia.obspm.fr¡/A_ , 2011, pp. 16
  * [15] Emiel H. Por et al.  “High Contrast Imaging for Python (HCIPy): an open-source adaptive optics and coronagraph simulator”  In _Adaptive Optics Systems VI_ 10703, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 2018, pp. 1070342  DOI: 10.1117/12.2314407
  * [16] Lisa A. Poyneer, Bruce A. Macintosh and Jean-Pierre Véran  “Fourier transform wavefront control with adaptive prediction of the atmosphere”  In _Journal of the Optical Society of America A_ 24.9, 2007, pp. 2645  DOI: 10.1364/JOSAA.24.002645
  * [17] Lisa A. Poyneer et al.  “Laboratory demonstration of the prediction of wind-blown turbulence by adaptive optics at 8 kHz with use of LQG control”  In _Applied Optics_ 62.8, 2023, pp. 1871  DOI: 10.1364/AO.474730
  * [18] Lisa A. Poyneer et al.  “Performance of the Gemini Planet Imager’s adaptive optics system”  In _Applied Optics_ 55.2, 2016, pp. 323  DOI: 10.1364/AO.55.000323
  * [19] Gaetano Sivo et al.  “First on-sky SCAO validation of full LQG control with vibration mitigation on the CANARY pathfinder”  In _Optics Express_ 22.19, 2014, pp. 23565  DOI: 10.1364/OE.22.023565
  * [20] Robin Swanson et al.  “Closed loop predictive control of adaptive optics systems with convolutional neural networks”  In _Monthly Notices of the Royal Astronomical Society_ 503.2, 2021, pp. 2944–2954  DOI: 10.1093/mnras/stab632
  * [21] Maaike van Kooten, Niek Doelman and Matthew Kenworthy  “Impact of time-variant turbulence behavior on prediction for adaptive optics systems”  In _Journal of the Optical Society of America A_ 36.5, 2019, pp. 731  DOI: 10.1364/JOSAA.36.000731
  * [22] Maaike A.. van Kooten et al.  “Predictive wavefront control on Keck II adaptive optics bench: on-sky coronagraphic results”  In _Journal of Astronomical Telescopes, Instruments, and Systems_ 8, 2022, pp. 029006  DOI: 10.1117/1.JATIS.8.2.029006
  * [23] Alison P. Wong et al.  “Predictive control for adaptive optics using neural networks”  In _Journal of Astronomical Telescopes, Instruments, and Systems_ 7, 2021, pp. 019001  DOI: 10.1117/1.JATIS.7.1.019001


# A Spatial Agent-Based Model for Preemptive Evacuation Decisions During
Typhoon

Rey C. Rodrigueza
Information and Communications Technology Department
Sorsogon State University - Bulan
Bulan, Sorsogon
<EMAIL_ADDRESS>
& Ma.Regina Justina E. Estuar
Ateneo Social Computing Science Laboratory
Ateneo de Manila University
Loyola Heights, Quezon City
<EMAIL_ADDRESS>

###### Abstract

Natural disasters continue to cause tremendous damage to human lives and
properties. The Philippines, due to its geographic location, is considered a
natural disaster-prone country experiencing an average of 20 tropical cyclones
annually. Understanding what factors significantly affect decision making
during crucial evacuation stages could help in making decisions on how to
prepare for disasters, how to act appropriately and strategically respond
during and after a calamity. In this work, an agent-based model for preemptive
evacuation decisions during typhoon is presented. In the model, civilians are
represented by households and their evacuation decisions were based from
calculated perceived risk. Also, rescuer and shelter manager agents were
included as facilitators during the preemptive evacuation process. National
and municipal census data were employed in the model, particularly for the
demographics of household agents. Further, geospatial data of a village in a
typhoon-susceptible municipality was used to represent the environment. The
decision to evacuate or not to evacuate depends on the agent’s perceived risk
which also depends on three decision factors: characteristics of the decision
maker (CDM); capacity related factors (CRF); and hazard related factors (HRF).
Finally, the number of households who decided to evacuate or opted to stay as
influenced by the model’s decision factors were determined during simulations.
Sensitivity analysis using linear regression shows that all parameters used in
the model are significant in the evacuation decision of household agents.

_Keywords_ Agent-Based Model $\cdot$ Typhoon $\cdot$ Preemptive Evacuation
Decisions

##  1 Introduction

Natural disasters are exacting an enormous damage with hundreds of thousands
of lives lost in the last decade alone (UNISDR, 2015). Asia appears to be the
mostly affected continent as in all disaster events, 44% of it were caused by
storms and floods which took 58% of the total deaths and 70% of the people
affected (Below and Wallemacq, 2018). According to World Risk Index 2018
(Heintze et al., 2018), the Philippines is the third country in the world with
the highest disaster risk. Almost ninety-nine percent of population in the
Philippines is potentially exposed to tropical typhoon, seventy-nine percent
of which to the most dangerous class of hazard (Pesaresi et al., 2017). This
alarming fact necessitates continuously finding solutions to reduce deaths and
damage during extreme typhoons and disasters given the fact that the
geographical location of the Philippines makes it vulnerable to natural
calamities or disasters. In doing so, loss of lives and properties could be
lessened, if not totally prevented.

Understanding how human behave and make decisions during disasters could help
in response and rescue efforts. Specifically, understanding factors used to
determine whether a unit evacuates or not may help in crafting policy and
corresponding implementation plans on how to prepare for disasters and how to
properly act and strategically respond during and after a calamity. Some
factors that have considerable effect to human behavior during disasters
should also be considered. Forewarning, for an instance, can do much to
mitigate consequence of disasters.

The Philippines’ weather state bureau PAGASA uses its Public Storm Warning
System to categorize typhoon intensity and warn the public (DOST-PAGASA, ).
The Public Storm Warning System has five Public Storm Warning Signals (PSWS)
as shown in Table 1. It also uses a set of color-coded rainfall advisories
comprising of 3 colors: yellow, orange, and red. Color yellow advisory warns
that flooding is expected in low lying areas and near river channels. Orange
advisory tells that flooding is threatening in low-lying areas and near river
channels. Finally, red rainfall advisory requires community response as severe
flooding is expected.

Table 1: Public Storm Warning System PSWS |  Lead Time (hrs) |  Winds (kph) |  Impacts of the Wind
---|---|---|---
#1 |  36 |  30 - 60 |  No damage to very light damage
#2 |  24 |  61 - 120 |  Light to moderate damage
#3 |  18 |  121 - 170 |  Moderate to heavy damage
#4 |  12 |  171 - 220 |  Heavy to very heavy damage
#5 |  12 |  $>$ 220 |  Very heavy to widespread damage

Typhoon causes several hazards like strong winds, flooding and landslides due
to incessant rain and storm surge prompting at-risk households to evacuate to
safe places. Most evacuations during typhoons are conducted in areas that are
prone to flooding. Household’s decision to evacuate is based on several
significant factors such as characteristics of the decision-maker in the
household, capacity-related factors and hazard-related factors (Lim et al.,
2016; Medina and Moraca, 2016; Cahyanto et al., 2014).

Identified decision maker’s attributes includes gender, age, level of
education, presence of young children, presence of elderly, presence of pet,
income level, house ownership, length of stay in residence and number of
household members (Lim et al., 2016; Cahyanto et al., 2014). Housing type and
previous typhoon experience on one hand are classified as capacity-related
factors (Whitehead et al., 2000; Hasan et al., 2010; Fu et al., 2006).
Finally, there are hazard-related factors that affects evacuation decision:
distance to typhoon, typhoon forward speed, typhoon wind speed, presence of
evacuation message, possibility of flooding, time of the day, source of notice
for evacuation and type of evacuation notice received (Hasan et al., 2010; Fu
et al., 2006).

Results on model estimation for partial and full household-level evacuation
decision model (Lim et al., 2016) for a large city in the Philippines indicate
that there is a higher probability of not evacuating if the head of the family
is male as compared to a female which is in consonance with the findings of
other studies (Horney et al., 2010; Lindell et al., 2005; Morrow and Gladwin,
2005); there is a higher probability that the households living in a house
with more than one floor level will not evacuate; there is more likelihood
that households will evacuate if the level of house damage is high; households
who own the house are more likely to evacuate compared to those who are just
renting; households with homes made from concrete materials have a higher
probability of staying at home compared to households whose houses are made of
wood; if the source of evacuation warning came from the authorities then the
households are more likely to evacuate compared to when the source of
information came from relatives or friends; in the presence of small children,
households are likely to evacuate which supports the results of other studies
(Cahyanto et al., 2014; Dash, 2002); and households that are located more than
10 meters away from the source of flood hazard have a high probability of
evacuating.

Agent-based modeling (ABM) is an approach to evaluate complex systems where
independent and interacting agents makes up its domain. It has also been used
to analyze systems that incorporate human behavior during decision-making
(Macal and North, 2005).

Many studies have used agent-based modeling for different purposes. It has
been applied to simulate disaster scenarios such as simulation of crowd
evacuation in a fire disaster (Wagner and Agrawal, 2014), flood incident
management (Dawson et al., 2011), and, human behavior during earthquake and
tsunami (D’Orazio et al., 2014; Mas et al., 2012) . However, very little study
has been gathered so far that used ABM to understand human behavior during
typhoon evacuation. This study attempts to model preemptive evacuation
decisions during typhoon.

##  2 Methodology

###  2.1 Overview of the area under study

Using ABM, a model was developed for San Vicente, a village in the
municipality of Bulan which is located at south-western most tip of the island
of Luzon, Philippines. The municipality of Bulan is almost always affected by
typhoons, experiencing an average of three to four typhoons per year. Hazard
maps of the municipality shows that it is prone to flooding, and San Vicente
is one of the most affected villages during typhoon (Mun, 2019).

###  2.2 ABM Development

The model environment is a GIS data of San Vicente village. The village is
mapped in OpenStreetMap.org and then the corresponding shapefiles are
downloaded from extract.bbbike.org. These shapefiles are collections of road
networks, areas, buildings, waterways, and points. The shapefiles are then
preprocessed using QGIS 3.4 before importing to GAMA Platform V1.8.0
(Taillandier et al., 2019) where the model is implemented. We opted to use
GAMA over other modeling platforms to build this model because integrating GIS
data in the model using GAMA is seamless.

####  2.2.1 Population Data

The population data came from two sources: 2015 census and housing data from
the Philippine Statistics Authority and local (Municipal) census data and
disaster report from Municipal Disaster Risk Reduction and Management Office
(MDRRMO) of Bulan.

####  2.2.2 Model Description

There are three types of human agents in the model. These are households,
rescuers, and shelter managers. The household agents are characterized by the
attributes of the head of the household: location, gender, income level, level
of education, has small kids in the household, has elderly in the household,
has household member with physical disability, type of house ownership, years
of residency, past typhoon experience, and perceived risk. Rescuer agents are
characterized by their location and capacity. On the other hand, a Shelter
Manager agent is characterized by the capacity of the evacuation shelter it is
assigned to. Typhoon wind intensity, rainfall severity, and time of day are
the exogenous factors/drivers of the model. Each household agent is assigned
to a building(house) and each evacuation center has an assigned shelter
manager agent. Rescuer agents are traced through their location coordinates.
Figure 1 shows the model during one of the simulation runs.

Figure 1: Model simulation showing household agents (white small circles)
evacuating.

####  2.2.3 Model Parameters and Initialization

Table 2 presents the model parameters with their type, definitions and initial
values. Other parameters of the model that are used for computing the values
of decision factors are presented in Tables 3, 4, and 5. The initialization
for household agents is always the same and the initial values for rescuers
and shelter managers are fixed. The initial values are based on data. These
data are sourced from national and local census and interviews from municipal
disaster managers.

Table 2: Model parameters and corresponding initial values Parameter |  Type |  Definition |  Initial Value
---|---|---|---
nb_households |  input (constant) |  number of household agents |  570
nb_rescuers |  input (constant) |  number of rescuer agents |  15
nb_sheltermanagers |  input (constant) |  number of shelter manager agents |  4
hhPerDistance |  input (constant) |  household agent perception distance in meters |  50
resPerDistance |  input (constant) |  rescuer agent perception distance in meters |  50
shelPerDistance |  input (constant) |  shelter manager agent perception distance in meters |  50
stormSeverity |  input (discrete) |  typhoon wind intensity |  1
rainfallSeverity |  input (discrete) |  rainfall intensity |  0.25
threshold |  input (discrete) |  threshold of evacuation decision |  0.7
timeOfDay |  input (discrete) |  time of day (daytime or night time) |  0.5
weight_CDM |  input (continuous) |  weight for Characteristics of Decision Maker (CDM) decision factor |  0.1
weight_HRF |  input (continuous) |  weight for Hazard Related Factors (HRF) decision factor |  0.1
weight_CRF |  input (continuous) |  weight for Capacity Related Factor (CRF) decision factor |  0.1
evacuated |  outcome |  number of household agents who evacuated |  0

At time t = 0, 570 household agents are initialized with their characteristics
values from a CSV file and then randomly assigned to houses/buildings. There
are also 15 rescuer agents randomly placed to their initial positions and 4
shelter managers assigned to four evacuation centers inside the village.
Rescuers and shelter manager’s parameters are initialized as well. Shapefiles
are imported to form the road networks, evacuation points, rescuer’s starting
points, buildings, as well as water bodies. The order in which household
agents decides to evacuate is also random as it depends on their individual
perceived risk value. The direction in which the rescuers are headed to inform
households about preemptive evacuation is also random.

Table 3: Characteristics of Decision Maker and Coding Format Head of household |  headHousehold |  Level of education |  educLevel |  Income level |  incomeLevel
---|---|---|---|---|---
Male |  0.5 |  College |  0.25 |  High |  0.25
Female |  1.0 |  High school |  0.5 |  Middle |  0.5
|  |  Grade school |  1.0 |  Low |  1.0
House ownership |  houseOwnership |  Presence of young children |  hasChildren |  Presence of elderly |  hasElderly
Owns the house |  0.5 |  No |  0.0 |  No |  0.0
Renting |  1.0 |  Yes |  1.0 |  Yes |  1.0
|  |  Presence of household member with disability |  withDisability |  Length of stay in the residence |  yearsOfResidency
|  |  No |  0.0 |  $>$ 10 years |  0.5
|  |  Yes |  1.0 |  $\leq$ 10 years |  1.0
Table 4: Hazard-Related Factors and and Coding Format Storm severity |  stormSeverity |  Proximity to the source of hazard |  proximityToHazard |  Time of day |  timeOfDay
---|---|---|---|---|---
PSWS 1 |  0.25 |  Far |  0.25 |  Daytime |  0.5
PSWS 2 |  0.5 |  Near |  0.5 |  Night time |  1.0
PSWS 3 |  1.0 |  Within |  1.0 |  |
Rainfall severity |  rainfallSeverity |  Source of evacuation warning |  sourceOfWarning |  |
Yellow |  0.25 |  Friends |  0.25 |  |
Orange |  0.5 |  Media |  0.5 |  |
Red |  1.0 |  Authorities |  1.0 |  |
Table 5: Capacity-Related Factors and and Coding Format House quality |  houseQuality |  House floor levels |  floorLevels |  Past typhoon experience |  typhoonExperience
---|---|---|---|---|---
Concrete |  0.25 |  $>$1 level |  0.5 |  Yes |  0.5
Wood |  0.5 |  1 level |  1.0 |  No |  1.0
Light materials |  1.0 |  |  |  |

####  2.2.4 Process Overview

Before the typhoon hits the village, rescuer agents roams around the village
and informs the household agents about the impending typhoon, then announces
preemptive evacuation. During this time, household agents computes their
perceived risk. If the computed perceived risk is higher than the percentage
of evacuation decision threshold to the highest possible perceived risk score
then the household agent evacuates to the nearest evacuation shelter inside
the village or outside of it. However, if the perceived risk value is lower
than the threshold value then the household agent stays (does not evacuate).
Shelter Manager agents accepts evacuees in evacuation centers and keeps check
on the shelter’s maximum capacity. If maximum capacity of the shelter is
reached then incoming evacuees will be directed to other evacuation shelters.

####  2.2.5 Individual Decision-Making

The household agents’ main objective is to survive the onslaught of typhoon.
The rescuer agents aim is to achieve zero casualty by informing household
agents about the typhoon in advance and announcing preemptive evacuation.
Shelter manager agents manages order and capacity of the evacuation center.
Household agents computes their perceived risk value in order to make a
decision. The perceived risk depends on three decision factors:
characteristics of decision maker; capacity-related factor; and hazard-related
factor. Bounded rationality is incorporated in the perceived-risk equation by
adding a random value from a known range of values. Household agents adapt to
exogenous state variables since these variables have direct effect on their
risk perception and individual objectives. Typhoon wind intensity, rainfall
severity, and time of day are the exogenous factors/drivers of the model. For
instance, an increasing public storm warning signal increases the hazard-
related factor value which is a component of the perceived risk equation for
household agents. In this model, spatial aspects play a role in the decision
process of the agents. Since the model environment is a GIS data, the distance
between objects is based from the real world. An agent that is near the source
of hazard such as river will have a higher risk perception compared to others
who are at a distance. Evacuating household agents also will choose the
nearest available evacuation center. Additionally, temporal aspects play a
role in the decision process as household agents perceives night time as a
more risky situation compared to daytime during typhoon.

Interactions between agents depend on spatial distances. Each agent type has a
perception area in a form of a circle based from perception radius which is
set at 50 meters. An agent can be sensed by another agent if it is inside the
perception area (which is in a form a circle) of the former. Figure 2 shows
the perception distances of the three agents in the model.

Figure 2: Perception Distances of Household, Rescuer and Shelter Manager
Agents.

###  2.3 Computation of the Perceived Risk

Perceived Risk (PR) of each household agents is based on three decision
factors: characteristics of the decision maker (CDM), hazard-related factor
(HRF), and capacity-related factor (CRF). PR is computed using equation 4.
Equation 1 computes the value for the characteristics of the decision maker
(CDF) in the household. The computation is performed one-time only during
simulation. HRF and CRF values are computed by equations 2 and 3 respectively.

$CDM\ =\ headHousehold\ +\ incomeLevel\ +\ educLevel+\ hasChildren\ +\
hasElderly\ \\\ +\ withDisability\ +\ houseOwnership\ +\ yearsOfResidency$ (1)
$HRF=stormSeverity\ +\ rainfallSeverity\ +\ proximityToHazard+\
sourceOfWarning\ \\\ +\ timeOfDay$ (2) $CRF=houseQuality\ +\ floorLevels\ +\
typhoonExperience$ (3) $Perceived\ Risk=(CDM*weight\\_CDM)\ +\
(HRF*weight\\_HRF)\ +\ (CRF*weight\\_CRF)\ +\ \varepsilon$ (4)

where $\varepsilon$ is the bounded rationality with value between 0.0 to 0.05.
The highest possible score for perceived risk is computed as

$Highest\ Possible\ Score=8\ *\ weight\\_CDM\ +\ 3*weight\\_CRF\ +\
5*weight\\_HRF$ (5)

To compute for the highest possible score of perceived risk, the maximum risk
value of the factors for each decision factors are added, then multiplied to
the assigned weight. For example, CDM decision factors is composed of eight
factors where the maximum risk value for each factor is 1.0. These maximum
risk values are then added which totalled to 8, then it is multiplied to the
assigned CDM weight. The same method is applied for HRF and CDF decision
factors. Finally, if household agent’s perceived risk value is higher than the
percentage of evacuation decision threshold to the highest possible perceived
risk score, then that household agent will evacuate.

###  2.4 Experiment Design

The model have seven experiment input parameters which includes stormSeverity,
rainfallSeverity, weight_CDM, weight_HRF, weight_CRF, and threshold. Five of
the seven parameters have constant increment step values while the other two.
rainfallSeverity and timeOfDay takes on the values 0.25, 0.5 , 1.0 and 0.5,
1.0 respectively. Table 6 show the minimum and maximum values of the input
parameters. There is only one outcome parameter for this model, the parameter
evacuated. This parameter counts the number of household agents who evacuated
upon assessing its perceived risk.

Table 6: Input Parameters Input Parameter |  min |  max |  step
---|---|---|---
stormSeverity |  1 |  2 |  1
rainfallSeverity |  0.25 |  1.0 |
timeOfDay |  0.5 |  1.0 |
threshold |  0.7 |  0.9 |  0.1
weight_CDM |  0.1 |  0.8 |  0.1
weight_HRF |  0.1 |  0.8 |  0.1
weight_CRF |  0.1 |  0.8 |  0.1

The model’s parameter space has 18,432 sets of possibilities using exhaustive
exploration method: stormSeverity (2) * rainfallSeverity (3) * timeOfDay (2) *
weight_CDM (8) * weight_HRF (8) * weight_CRF (8) * threshold (3) = 18,432.
Each possibility is replicated 10 times thus generating 184,320 total
simulations for the batch experiment. It is to be noted, however, that for a
certain parameter combinations or possibility to be considered for simulation
run, added weights of CDM, HRF and CRF must total to 1.0. If the total weight
is less than 1.0 then the current possibility simulation is skipped and
proceeds to the next one. Of the 18,432 sets of possibilities, only 6,480 are
considered and the rest are skipped.

###  2.5 Sensitivity Analysis

Sensitivity analysis is used to quantify which model inputs are more
influential for the model output. To determine significant predictors of
evacuation decision, we used the linear regression function lm() of the R
language on the output data from the batch experiment. The independent
variables are the seven sampled input parameters while the dependent variable
is the sole model outcome, the number of household agents who evacuated.

##  3 Results and Discussion

Sensitivity analysis shows that all parameters used in the model (weight_CDM,
weight_HRF, weight_CRF, stormSeverity, rainfallSeverity, timeOfDay, threshold)
are significant in the evacuation decision of household agents having
significance levels of $<2e^{-16}$. The $p$-value is also $<2e^{-16}$.
Findings also show that the total evacuation decision is more sensitive to
threshold and weights assigned in capacity related factors, specifically type
of house, floor level and past typhoon experience.

(a) two,orange,nighttime

(b) two,orange,nighttime

(c) two,orange,nighttime

Figure 3: Simulation at (a) 70%, (b) 80%, and (c) 90% threshold of evacuation
decision

Figure 3 shows selected three simulation results from the batch experiment
conducted. The label for each graph follows the format {storm severity (public
storm warning signal), rainfall severity, time of day}. The three colored
graph lines represents values for weight_CRF (red), weight_HRF (green), and
weight_CDM_ (blue). The result of one of the three simulation runs, simulation
(c), is near to the real-world data based on preemptive evacuation happened
during a typhoon that passed in San Vicente village on 2019 where 66
households preemptively evacuated. Simulation (c) result is near to the real-
world data based on preemptive evacuation happened during a typhoon that
passed in San Vicente village on 2019 where 66 households preemptively
evacuated. This indicates that at this typhoon condition, some household
decision makers might have reached their evacuation decision threshold and
finally decided that the best action is to evacuate. Not doing so would be
catastrophic at their end. It is also noteworthy that for all threshold
settings, an increase in the weight for capacity-related factors also
increases the number of household agents who decided to evacuate. On the other
hand, the effect of increasing the weights for characteristics of decision
makers and hazard-related factors on the number of evacuees is not evident.

The batch experiment can actually be segmented into three parts. The first
part is where the threshold of evacuation decision is at 70% and the other two
parts are for 80% and 90% thresholds respectively. At 70% threshold, the
number of evacuees is high and even reaches to 570 which is the total number
of household agents even when the actual risk is low, such as conditions
having PSWS 1, yellow rainfall, and day time. The number of evacuees gradually
decreased when the threshold is at 80%, but the numbers are still high
compared to real-world data. However, at 90% threshold, evacuees only started
to appear and increase at condition {two,red,daytime}. The number of evacuees
is at nearest to real-world data at condition {two,orange,nighttime},
indicating that at this typhoon condition, some household decision makers
might have reached their evacuation decision threshold and finally decided
that the best action to do is to evacuate. Not doing so would be catastrophic
at their end. It is also important to note that for all threshold settings, an
increase in the weight for capacity-related factors also increases the number
of household agents who decided to evacuate. On the other hand, the effect of
increasing the weights for characteristics of decision makers and hazard-
related factors on the number of evacuees is not that evident. Due to limited
space, only the simulation result at 80% threshold is presented in this paper
and is shown in Figure 4.

(a) one,yellow,daytime

(b) one,orange,daytime

(c) one,red,daytime

(d) one,yellow,nighttime

(e) one,orange,nighttime

(f) one,red,nighttime

(g) two,yellow,daytime

(h) two,orange,daytime

(i) two,red,daytime

(j) two,yellow,nighttime

(k) two,orange,nighttime

(l) two,red,nighttime

Figure 4: Simulation at 80% threshold of evacuation decision.

##  4 Conclusion

This paper presented a spatial agent-based model for investigating preemptive
evacuation decisions during typhoon. GIS data of a typhoon-susceptible village
in south-western part of the Philippines is used as the environment of the
model and the data utilized is based on national and local census. Among the
three human agents in the model, the focus is more on the heterogeneous
household agents and their perceived risk during typhoon. The perceived risk
of household agents which depends on three decision factors is the major
determinant whether a household agent preemptively evacuates or not. Model
simulation was performed using exhaustive exploration method on the set of
input parameters. A single output variable determines the number of evacuation
decisions.

The behavior of the agents was reviewed. It was seen that among the perceived
risk decision factors, capacity-related factors contributes significantly in
the decision of household agents to evacuate. The characteristics of the
decision maker in a household and hazard-related factors also contributes but
not as significant as capacity-related factors. Another deciding element
observed is the threshold of evacuation decision. This threshold value is the
assumed perceived limit of what a household decision maker can hold before
finally deciding that the best action to do is to evacuate. These findings
could provide insights to disaster managers and health emergency officers as
to why households decides to preemptively evacuate while others opt to stay
and wait for weather developments. It could also give initial understanding on
how perceived risk and its decision factors directly affects how people behave
during disaster such as typhoon.

## Acknowledgement

The authors would like to thank the Municipal Disaster Risk Reduction and
Management Office (MDRRMO) of LGU Bulan, Sorsogon for providing the necessary
information used in this study and the Ateneo Social Computing Science
Laboratory, Ateneo de Manila University for providing resources for the
conduct of this study.

## References

  * UNISDR [2015] UNISDR.  Global assessment report on disaster risk reduction making development sustainable: the future of disaster risk management.  Technical report, Geneva, Switzerland: UNISDR, 2015.
  * Below and Wallemacq [2018] R. Below and P. Wallemacq.  Annual disaster statistical review 2017.  Technical report, Centre for Research on the Epidemiology of Disasters (CRED)-Université Catholique de Louvain, Brussels, Belgium, 2018.
  * Heintze et al. [2018] Hans-Joachim Heintze, Lotte Kirch, Barbara Kuppers, Holger Mann, Frank Mischo, Peter Mucke, Tanja Pazdzierny, Ruben Prutz, Katrin Radtke, Friederike Strube, and Daniel Weller.  Worldriskreport 2018.  Technical report, Bundnis Entwicklung Hilft, 2018.
  * Pesaresi et al. [2017] Martino Pesaresi, Daniele Ehrlich, Thomas Kemper, Alice Siragusa, Aneta J Florczyk, Sergio Freire, and Christina Corbane.  Atlas of the human planet 2017. global exposure to natural hazards.  Technical report, Publications Office of the European Union, 2017.
  * [5] DOST-PAGASA.  Revised public storm warning system.  URL https://bagong.pagasa.dost.gov.ph/learning-tools/public-storm-warning-signal.
  * Lim et al. [2016] Ma Bernadeth B Lim, Hector R Lim, Mongkut Piantanakulchai, and Francis Aldrine Uy.  A household-level flood evacuation decision model in quezon city, philippines.  _Natural Hazards_ , 80(3):1539–1561, 2016.
  * Medina and Moraca [2016] Michael Arieh P Medina and Jackielyn M Moraca.  Should i stay or should i go? determinants of evacuation upon flood warning among households in a flood prone area in bukidnon, philippines.  _International Letters of Natural Sciences_ , 50, 2016.
  * Cahyanto et al. [2014] Ignatius Cahyanto, Lori Pennington-Gray, Brijesh Thapa, Siva Srinivasan, Jorge Villegas, Corene Matyas, and Spiro Kiousis.  An empirical evaluation of the determinants of tourist’s hurricane evacuation decision making.  _Journal of Destination Marketing & Management_, 2(4):253–265, 2014.
  * Whitehead et al. [2000] John C Whitehead, Bob Edwards, Marieke Van Willigen, John R Maiolo, Kenneth Wilson, and Kevin T Smith.  Heading for higher ground: factors affecting real and hypothetical hurricane evacuation behavior.  _Global Environmental Change Part B: Environmental Hazards_ , 2(4):133–142, 2000.
  * Hasan et al. [2010] Samiul Hasan, Satish Ukkusuri, Hugh Gladwin, and Pamela Murray-Tuite.  Behavioral model to understand household-level hurricane evacuation decision making.  _Journal of Transportation Engineering_ , 137(5):341–348, 2010.
  * Fu et al. [2006] Haoqiang Fu, Chester Wilmot, and Earl Baker.  Sequential logit dynamic travel demand model and its transferability.  _Transportation Research Record: Journal of the Transportation Research Board_ , (1977):17–26, 2006.
  * Horney et al. [2010] Jennifer A Horney, Pia DM MacDonald, Marieke Van Willigen, Philip R Berke, and Jay S Kaufman.  Individual actual or perceived property flood risk: Did it predict evacuation from hurricane isabel in north carolina, 2003?  _Risk Analysis: An International Journal_ , 30(3):501–511, 2010.
  * Lindell et al. [2005] Michael K Lindell, Jing-Chein Lu, and Carla S Prater.  Household decision making and evacuation in response to hurricane lili.  _Natural Hazards Review_ , 6(4):171–179, 2005\.
  * Morrow and Gladwin [2005] Betty H. Morrow and Huge Gladwin.  Hurricane ivan behavioral analysis, 2004 hurricane assessments.  Technical report, Federal Emergency Management Agency and U.S. Army Corps of Engineers, Washington, DC, 2005.
  * Dash [2002] Nicole Dash.  Decision-making under extreme uncertainty: Rethinking hazard related perceptions and action.  2002\.
  * Macal and North [2005] Charles M Macal and Michael J North.  Tutorial on agent-based modeling and simulation.  In _Simulation conference, 2005 proceedings of the winter_ , pages 14–pp. IEEE, 2005.
  * Wagner and Agrawal [2014] Neal Wagner and Vikas Agrawal.  An agent-based simulation system for concert venue crowd evacuation modeling in the presence of a fire disaster.  _Expert Systems with Applications_ , 41(6):2807–2815, 2014.
  * Dawson et al. [2011] Richard J Dawson, Roger Peppe, and Miao Wang.  An agent-based model for risk-based flood incident management.  _Natural hazards_ , 59(1):167–189, 2011.
  * D’Orazio et al. [2014] Marco D’Orazio, Luca Spalazzi, Enrico Quagliarini, and Gabriele Bernardini.  Agent-based model for earthquake pedestrians’ evacuation in urban outdoor scenarios: Behavioural patterns definition and evacuation paths choice.  _Safety science_ , 62:450–465, 2014.
  * Mas et al. [2012] Erick Mas, Anawat Suppasri, Fumihiko Imamura, and Shunichi Koshimura.  Agent-based simulation of the 2011 great east japan earthquake/tsunami evacuation: An integrated model of tsunami inundation and evacuation.  _Journal of Natural Disaster Science_ , 34(1):41–57, 2012.
  * Mun [2019] _Contingency Plan for Hydro-meteorological Hazards and its Secondary Hazards of Flood, Storm Surge and Landslides: A Dynamic Ten-Year Plan (2019-2028)_.  Municipal Disaster Risk Reduction and Management Council of LGU Bulan, Philippines, 2019.
  * Taillandier et al. [2019] Patrick Taillandier, Benoit Gaudou, Arnaud Grignard, Quang-Nghi Huynh, Nicolas Marilleau, Philippe Caillou, Damien Philippon, and Alexis Drogoul.  Building, composing and experimenting complex spatial models with the gama platform.  _GeoInformatica_ , 23(2):299–322, 2019.