* [39] Danish Health Authorities  “Vejledning og Instruktion Vedrørende Sundhedsplejerskeordningen [Guidance and Instructions on the Infant Health Nurse Visiting Program]”, 1954
  * [40] Danish Health Authorities  “Vejledning og Instruktion Vedrørende Sundhedsplejerseordningen [Guidance and Instructions on the Infant Health Nurse Visiting Program]”, 1961  URL: https://www.xn--betnkninger-c9a.dk/wp-content/uploads/2021/02/573.pdf
  * [41] Det Statistiske Departement  “Spædbørnsdødelighed 1931-1960 [Infant Mortality 1931-1960]”  Det Statistiske Departement, 1964
  * [42] Alan Dolan and Joan Aldous  “Networks and Algorithms: An Introductory Approach”  Chichester: John Wiley & Sons Ltd., 1993
  * [43] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit and Neil Houlsby  “An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale”  In _ICLR_ , 2021
  * [44] Miroslav Dudík, John Langford and Lihong Li  “Doubly robust policy evaluation and learning”  In _Proceedings of the 28th International Conference on International Conference on Machine Learning_ , 2011, pp. 1097–1104
  * [45] Valentina Duque, Maria Rosales-Rueda and Fabio Sanchez  “How do early-life shocks interact with subsequent human-capital investments? Evidence from administrative data”  In _IZA World of Labor Conference_ , 2018
  * [46] Chloe N. East, Sarah Miller, Marianne Page and Laura R. Wherry  “Multi-generational impacts of childhood access to the safety net: Early life exposure to Medicaid and the next generation’s health”  In _American Economic Review_ 113.1, 2023, pp. 98–135
  * [47] Martin Ester, H.. Kriegel, Jorg Sander and Xiaowei Xu  “A density-based algorithm for discovering clusters in large spatial databases with noise”  In _Proc. 2nd Int. Conf. Knowledge Discovery and Data Mining (KDD’96)_ AAAI Press, 1996, pp. 226–231
  * [48] Vivek F Farias and Andrew A Li  “Learning preferences with side information”  In _Management Science_ 65.7 INFORMS, 2019, pp. 3131–3149
  * [49] Edgar C Fieller, Herman O Hartley and Egon S Pearson  “Tests for rank correlation coefficients. I”  In _Biometrika_ 44.3/4 JSTOR, 1957, pp. 470–481
  * [50] Lester Randolph Ford  “Network flow theory”  Rand Corporation Santa Monica, CA, 1956
  * [51] Lester Randolph Ford Jr and Delbert Ray Fulkerson  “Flows in Networks”  N.J.: Princeton University Press, 1962
  * [52] Anders Forsdahl  “Are poor living conditions in childhood and adolescence an important risk factor for arteriosclerotic heart disease?”  In _British Journal of Preventive and Social Medicine_ 31.2, 1979, pp. 91–95
  * [53] Yoshimi Fukuoka, Mo Zhou, Eric Vittinghoff, William Haskell, Ken Goldberg and Anil Aswani  “Objectively Measured Baseline Physical Activity Patterns in Women in the mPED Trial: Cluster Analysis”  In _JMIR Public Health and Surveillance_ 4.1 JMIR Publications Inc., Toronto, Canada, 2018, pp. e9138
  * [54] Wei Gao and Russ Tedrake  “FilterReg: Robust and Efficient Probabilistic Point-Set Registration Using Gaussian Filter and Twist Parameterization”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019, pp. 11095–11104
  * [55] Teevrat Garg, Maulik Jagnani and Vis Taraz  “Temperature and human capital in India”  In _Journal of the Association of Environmental and Resource Economists_ 7.6 The University of Chicago Press Chicago, IL, 2020, pp. 1113–1150
  * [56] R. Geetha, T. Thilagam and T. Padmavathy  “Effective offline handwritten text recognition model based on a sequence-to-sequence approach with CNN–RNN networks”  In _Neural Computing and Applications_ 33.17, 2021, pp. 10923–10934
  * [57] Andrew V Goldberg  “A new max-flow algorithm”, 1985
  * [58] Andrew V Goldberg  “Efficient graph algorithms for sequential and parallel computers”, 1987
  * [59] Andrew V Goldberg and Robert E Tarjan  “A new approach to the maximum-flow problem”  In _Journal of the ACM (JACM)_ 35.4 ACM New York, NY, USA, 1988, pp. 921–940
  * [60] Andrew V Goldberg and Robert E Tarjan  “Finding minimum-cost circulations by canceling negative cycles”  In _Journal of the ACM (JACM)_ 36.4 ACM New York, NY, USA, 1989, pp. 873–886
  * [61] Ian J. Goodfellow, Yaroslav Bulatov, Julian Ibarz, Sacha Arnoud and Vinay Shet  “Multi-digit Number Recognition from Street View Imagery using Deep Convolutional Neural Networks”  In _arXiv:1312.6082_ , 2013
  * [62] Andrew Goodman-Bacon  “Public insurance and mortality: Evidence from Medicaid implementation”  In _Journal of Political Economy_ 126.1 University of Chicago Press Chicago, IL, 2018, pp. 216–262
  * [63] Margaret A Hamburg and Francis S Collins  “The Path to Personalized Medicine”  In _New England Journal of Medicine_ 363.4 Mass Medical Soc, 2010, pp. 301–304
  * [64] James J. Heckman, Seong Hyeok Moon, Rodrigo Pinto, Peter A. Savelyev and Adam Yavitz  “The rate of return to the HighScope Perry Preschool Program”  In _Journal of Public Economics_ 94.1-2 Elsevier, 2010, pp. 114–128
  * [65] Keisuke Hirano and Jack R Porter  “Asymptotics for statistical treatment rules”  In _Econometrica_ 77.5 Wiley Online Library, 2009, pp. 1683–1701
  * [66] Keisuke Hirano and Jack R Porter  “Asymptotic analysis of statistical decision rules in econometrics”  In _Handbook of econometrics_ 7 Elsevier, 2020, pp. 283–354
  * [67] Jonas Hjort, Mikkel Sølvsten and Miriam Wüst  “Universal Investment in Infants and Long-Run Health: Evidence from Denmark’s 1937 Home Visiting Program”  In _American Economic Journal: Applied Economics_ 9.4, 2017, pp. 78–104  DOI: 10.1257/app.20150087
  * [68] Lauren Hoehn-Velasco  “The Long-term Impact of Preventative Public Health Programs”  In _The Economic Journal_ 131.634 Oxford University Press, 2021, pp. 797–826
  * [69] Hilary Hoynes, Diane Whitmore Schanzenbach and Douglas Almond  “Long-run impacts of childhood access to the safety net”  In _American Economic Review_ 106.4, 2016, pp. 903–34
  * [70] Den Danske Almindelige Lægeforenings Hygiejnekomite J.  “Vejledning vedrørende gennemførelsen af lov nr. 85 af 31. marts 1937 om bekæmpelse af sygelighed og dødelighed blandt børn i det første leveår [Guidance on the implementation of the act no. 85 of the 31st March, 1937 on Combating Morbidity and Mortality among Children in the First Year of Life]”  In _Ugeskrift for Laeger_ 100.10, 1938
  * [71] Kewal K Jain  “Personalized Medicine”  In _Current Opinion in Molecular Therapeutics_ 4.6, 2002, pp. 548–558
  * [72] Torben Skov Dyg Johansen  “timm-sequence-net” Python library to be made open-source.  In _GitHub repository_ GitHub, https://github.com/TorbenSDJohansen/timm-sequence-net, 2023
  * [73] Nathan Kallus  “Balanced Policy Evaluation and Learning”  In _Advances in Neural Information Processing Systems_ 31 Curran Associates, Inc., 2018
  * [74] Nathan Kallus and Angela Zhou  “Confounding-Robust Policy Improvement”  In _Advances in Neural Information Processing Systems_ 31 Curran Associates, Inc., 2018
  * [75] Lei Kang, Pau Riba, Marçal Rusiñol, Alicia Fornés and Mauricio Villegas  “Pay Attention to What You Read: Non-recurrent Handwritten Text-Line Recognition”  In _Pattern Recognition_ 129, 2022, pp. 108766
  * [76] Alexander V Karzanov  “Determining the maximal flow in a network by the method of preflows”  In _Soviet Math. Dokl._ 15, 1974, pp. 434–437
  * [77] Edward S Kim, Roy S Herbst, Ignacio I Wistuba, J Jack Lee, George R Blumenschein Jr, Anne Tsao, David J Stewart, Marshall E Hicks, Jeremy Erasmus Jr and Sanjay Gupta  “The BATTLE trial: Personalizing therapy for lung cancer”  In _Cancer Discovery_ 1.1 AACR, 2011, pp. 44–53
  * [78] Toru Kitagawa and Aleksey Tetenov  “Who should be treated? empirical welfare maximization methods for treatment choice”  In _Econometrica_ 86.2 Wiley Online Library, 2018, pp. 591–616
  * [79] Morton Klein  “A primal method for minimal cost flows with applications to the assignment and transportation problems”  In _Management Science_ 14.3 INFORMS, 1967, pp. 205–220
  * [80] William H Kruskal and W Allen Wallis  “Use of ranks in one-criterion variance analysis”  In _Journal of the American statistical Association_ 47.260 Taylor & Francis, 1952, pp. 583–621
  * [81] C. Lee and Simon Osindero  “Recursive recurrent nets with attention modeling for OCR in the wild”  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 2231–2239
  * [82] Charles F Manski  “Statistical treatment rules for heterogeneous populations”  In _Econometrica_ 72.4 Wiley Online Library, 2004, pp. 1221–1246
  * [83] Sarah Miller and Laura R. Wherry  “The long-term effects of early life Medicaid coverage”  In _Journal of Human Resources_ 54.3 University of Wisconsin Press, 2019, pp. 785–824
  * [84] Edward F Moore  “The shortest path through a maze”  In _Proc. of the International Symposium on the Theory of Switching_ , 1959, pp. 285–292  Harvard University Press
  * [85] Katta G Murty  “Network programming”  Englewood Cliffs, N.J.: Prentice-Hall, Inc., 1992
  * [86] Xinkun Nie and Stefan Wager  “Quasi-Oracle Estimation of Heterogeneous Treatment Effects”  In _Biometrika_ 108.2 Oxford University Press, 2021, pp. 299–319
  * [87] Hamid Noghanibehambari  “Intergenerational health effects of Medicaid”  In _Economics & Human Biology_ 45, 2022, pp. 101114  DOI: https://doi.org/10.1016/j.ehb.2022.101114
  * [88] David L. Olds, Charles R. Henderson, Robert Chamberlin and Robert Tatelbaum  “Preventing Child Abuse and Neglect: A Randomized Trial of Nurse Home Visitation”  In _Pediatrics_ 78.1 American Academy of Pediatrics, 1986, pp. 65–78  URL: http://pediatrics.aappublications.org/content/78/1/65
  * [89] David L. Olds, Charles R. Henderson, Robert Cole, John Eckenrode, Dennis Kitzmann, Lisa Pettitt, Kimberly Sidora, Pamela Morris and Jane Powers  “Long-term effects of nurse home visitation on childrens criminal and antisocial behavior: 15-year follow-up of a randomized controlled trial”  In _JAMA_ 280.14, 1998, pp. 1238–1244  DOI: 10.1001/jama.280.14.1238
  * [90] David L. Olds, Harriet Kitzman, Elizabeth Anson, Joyce A. Smith, Michael D. Knudtson, Ted Miller, Robert Cole, Christian Hopfer and Gabriella Conti  “Prenatal and infancy nurse home visiting effects on mothers: 18-year follow-up of a randomized trial”  In _Pediatrics_ 144 American Academy of Pediatrics, 2019, pp. e20183889
  * [91] Elissa M Ozanne, Rebecca Howe, Zehra Omer and Laura J Esserman  “Development of a personalized decision aid for breast cancer risk reduction and management”  In _BMC Medical Informatics and Decision Making_ 14.1 BioMed Central, 2014, pp. 1–8
  * [92] Maya Rossin-Slater and Miriam Wüst  “What Is the Added Value of Preschool for Poor Children? Long-Term and Intergenerational Impacts and Interactions with an Infant Health Intervention”  In _American Economic Journal: Applied Economics_ 12.3, 2020, pp. 255–286
  * [93] Nicholas J Schork  “Personalized medicine: Time for one-person trials”  In _Nature_ 520, 2015, pp. 609–611
  * [94] Zejiang Shen, Ruochen Zhang, Melissa Dell, Benjamin Charles Germain Lee, Jacob Carlson and Weining Li  “LayoutParser: A Unified Toolkit for Deep Learning Based Document Image Analysis”  In _International Conference on Document Analysis and Recognition_ Cham.: Springer, 2021, pp. 131–146
  * [95] Richard S Sutton and Andrew G Barto  “Reinforcement learning: An introduction”  MIT press, 2018
  * [96] Adith Swaminathan and Thorsten Joachims  “Batch learning from logged bandit feedback through counterfactual risk minimization”  In _The Journal of Machine Learning Research_ 16.1 JMLR. org, 2015, pp. 1731–1755
  * [97] Adith Swaminathan, Akshay Krishnamurthy, Alekh Agarwal, Miro Dudik, John Langford, Damien Jose and Imed Zitouni  “Off-policy evaluation for slate recommendation”  In _Advances in Neural Information Processing Systems_ 30, 2017
  * [98] Mingxing Tan and Quoc Le  “EfficientNetV2: Smaller Models and Faster Training”  In _Proceedings of the 38th International Conference on Machine Learning_ 139, Proceedings of Machine Learning Research  PMLR, 2021, pp. 10096–10106  URL: https://proceedings.mlr.press/v139/tan21a.html
  * [99] Stefan Wager and Susan Athey  “Estimation and inference of heterogeneous treatment effects using random forests”  In _Journal of the American Statistical Association_ 113.523 Taylor & Francis, 2018, pp. 1228–1242
  * [100] Laura R Wherry, Sarah Miller, Robert Kaestner and Bruce D Meyer  “Childhood medicaid coverage and later-life health care utilization”  In _Review of Economics and Statistics_ 100.2 MIT Press, 2018, pp. 287–302
  * [101] Ross Wightman  “PyTorch Image Models”  In _GitHub repository_ GitHub, https://github.com/rwightman/pytorch-image-models, 2019  DOI: 10.5281/zenodo.4414861
  * [102] Miriam Wüst  “Early interventions and infant health: Evidence from the Danish home visiting program”  In _Labour Economics_ 19.4 Elsevier, 2012, pp. 484–495
  * [103] Baqun Zhang, Anastasios A Tsiatis, Marie Davidian, Min Zhang and Eric Laber  “Estimating optimal treatment regimes from a classification perspective”  In _Statistics_ 1.1 Wiley Online Library, 2012, pp. 103–114
  * [104] Ying-Qi Zhao, Donglin Zeng, Eric B Laber, Rui Song, Ming Yuan and Michael Rene Kosorok  “Doubly robust learning for estimating individualized treatment with censored data”  In _Biometrika_ 102.1 Oxford University Press, 2015, pp. 151–168
  * [105] Yingqi Zhao, Donglin Zeng, A John Rush and Michael R Kosorok  “Estimating individualized treatment rules using outcome weighted learning”  In _Journal of the American Statistical Association_ 107.499 Taylor & Francis, 2012, pp. 1106–1118
  * [106] Mo Zhou, Yoshimi Fukuoka, Yonatan Mintz, Ken Goldberg, Philip Kaminsky, Elena Flowers and Anil Aswani  “Evaluating machine learning–based automated personalized daily step goals delivered through a mobile phone app: Randomized controlled trial”  In _JMIR mHealth and uHealth_ 6.1 JMIR Publications Inc., Toronto, Canada, 2018, pp. e9117
  * [107] Xin Zhou, Nicole Mayer-Hamblett, Umer Khan and Michael R Kosorok  “Residual weighted learning for estimating individualized treatment rules”  In _Journal of the American Statistical Association_ 112.517 Taylor & Francis, 2017, pp. 169–187
  * [108] Zhengyuan Zhou, Susan Athey and Stefan Wager  “Offline multi-action policy learning: Generalization and optimization”  In _Operations Research_ 71.1 INFORMS, 2023, pp. 148–183

## Appendices

###  A Extended Proof

This section provides additional details on the proofs of Theorems 2.3 and
2.5.

###### Extended Proof of Theorem 2.3.

For each $v_{1}\in V_{1}$ (treatment) and $v_{2}\in V_{2}$ (recipient), let
$D_{v_{2}}^{v_{1}}=x_{v_{1}v_{2}}$ with $x$ a feasible flow in $\mathcal{N}$.
To prove that $x$ corresponds to a valid solution of the optimal allocation
problem, it is sufficient to prove that any feasible flow $x$ respects the
constraints of the optimal allocation problem’s binary program representation.

Recall that a flow $x$ is feasible if, for all $uw\in A$, $l_{uw}\leq
x_{uw}\leq u_{uw}$ and, for all $v\in V$, $b(v)=b_{x}(v)$. Exactly one arc
enters each $v_{1}\in V_{1}$, namely $sv_{1}$, and its capacity is
$u_{sv_{1}}=m$. Since $b(v_{1})=b_{x}(v_{1})$, this implies that at most $m$
units of flow leaves $v_{1}$, and since $x$ is an integer flow this implies
that for each $v_{1}\in V_{1}$, at most $m$ of the arcs $v_{1}v_{2}$,
$v_{2}\in V_{2}$, have non-zero flow. This proves that a feasible flow
allocates any treatment to at most $m$ recipients.

Exactly one arc leaves each $v_{2}\in V_{2}$, namely $v_{2}t$, and its lower
bound and capacity are $l_{v_{1}t}=1$ and $u_{v_{2}t}=1$, respectively. Since
$b(v_{2})=b_{x}(v_{2})$ and $x$ is an integer flow, this implies that for each
$v_{2}\in V_{2}$, exactly one of the arcs $v_{1}v_{2}$, $v_{1}\in V_{1}$, have
non-zero flow. This proves that a feasible flow allocates exactly one
treatment to any recipient.

To prove the relationship between the cost of a feasible flow and the value of
a solution to the optimal allocation problem’s binary program representation,
note that the cost of any feasible flow $x$ is:

$\displaystyle\sum_{uw\in A}c_{uw}x_{uw}$ $\displaystyle=\sum_{v_{1}\in
V_{1}}\sum_{v_{2}\in V_{2}}c_{v_{1}v_{2}}x_{v_{1}v_{2}}$ Since $c_{uw}$ is 0
for all other arcs $\displaystyle=\sum_{v_{2}\in V_{2}}\sum_{v_{1}\in
V_{1}}c_{v_{1}v_{2}}x_{v_{1}v_{2}}$ Switching order of summations
$\displaystyle=\sum_{v_{2}\in V_{2}}\sum_{v_{1}\in
V_{1}}c_{v_{1}v_{2}}D_{v_{2}}^{v_{1}}$ Since
$x_{v_{1}v_{2}}=D_{v_{2}}^{v_{1}}$ for all $v_{1}\in V_{1}$, $v_{2}\in V_{2}$
$\displaystyle=-\sum_{v_{2}\in V_{2}}\sum_{v_{1}\in
V_{1}}Y_{v_{2}}^{v_{1}}D_{v_{2}}^{v_{1}}$ Using
$c_{v_{1}v_{2}}=-Y_{v_{2}}^{v_{1}}$

∎

To prove Theorem 2.5, we shall make use of the following algorithm [79]:

Algorithm 1 Cycle Cancelling Algorithm

A network $\mathcal{N}=(V,A,l,u,b,c)$.

A minimum cost feasible flow $x$ in $\mathcal{N}$.

procedure CycleCanceling($\mathcal{N}$)

Find a feasible flow $x$ in $\mathcal{N}$.

repeat

Search for a negative cycle $C$ in $\mathcal{N}(x)$.

if such a cycle $C$ is found then

Augment $x$ by $\delta(X)$ units along $C$.

end if

until no negative cycle is found

return $x$.

end procedure

###### Extended Proof of Theorem 2.5.

Let $n=|V|$ denote the number of vertices and $m=|A|$ the number of arcs in
the underlying digraph $D$ of $\mathcal{N}$. Using that the Bellman-Ford-Moore
algorithm allows us to check for the existence of a negative cycle in an
arbitrary digraph in time $\mathcal{O}(nm)$ [23, 50, 84], we may use the cycle
cancelling algorithm [79] to find a minimum cost feasible flow in a network
(Algorithm 1).656565This is due to the result that $x$ is a minimum cost
feasible flow in $\mathcal{N}$ if and only if $\mathcal{N}(x)$ contains no
directed cycle of negative cost [14]. However, in the absence of guidance as
to how to find the cycles to augment along, this algorithm is guaranteed to
work only when lower bounds, capacities, costs, and balance vectors are
integers, and may be exponential in the size of the input (i.e., not strongly
polynomial).

The problem with the “naïve” cycle cancelling algorithm is in its searching
step for a negative cycle. Without searching for an appropriate negative
cycle, the algorithm is not useful for the type of network we consider. To
confront this issue, we use the approach pioneered in [57] and generalized in
[59] of using preflows [76] for solving maximum flow problems. Using
generalizations due to [58, 60], this allows us to modify the searching step
of the algorithm in a way to obtain a strongly polynomial algorithm by always
augmenting along a cycle of minimum mean cost.666666The mean cost of a cycle
$C$ is defined as $c(C)/|A(C)|$. Specifically, we can find a minimum mean cost
cycle in time $\mathcal{O}(\log n)$, resulting in total time
$\mathcal{O}(nm^{2}(\log n)^{2})$.

Recall that we require any solution to be an integer flow.676767Otherwise some
recipient could be allocated different “shares” of different treatments, which
is incompatible with the binary program structure of the optimal allocation
problem. Here, we may use a convenient integrality property of minimum cost
flows, namely that given all integer lower bounds, capacities, and balance
vectors, there exists an integer minimum cost flow. To see why, we may start
by assuming, without loss of generality, that our initial feasible flow is an
integer flow.686868Either by using the initial allocation as the flow or by
exploiting the layered structure of $\mathcal{N}$ to easily find a feasible
integer flow in time $\mathcal{O}(m)$. At each iteration, we augment our flow
by $\delta(C)$, by which we mean the minimum residual capacity of any arc on
$C$ in $\mathcal{N}$, and thus our new flow after any iteration is
$x^{t+1}=x^{t}\oplus\delta(C)$. Since the residual capacity on any arc is
always an integer, $x^{t+1}$ is an integer flow so long as $x^{t}$ is, and
thus by induction we have an integer minimum cost flow.

Using the structure of the network, we can obtain the complexity in terms of
the inputs to the optimal allocation problem. Noting that $n=2+n_{1}+n_{2}$
and $m=n_{1}+n_{2}+n_{1}n_{2}$, we have:

$\displaystyle\mathcal{O}\left[nm^{2}(\log n)^{2}\right]$
$\displaystyle=\mathcal{O}\left[(2+n_{1}+n_{2})(n_{1}+n_{2}+n_{1}n_{2})^{2}\log^{2}(2+n_{1}+n_{2})\right]$
$\displaystyle=\mathcal{O}\left[(n_{1}+n_{2})(n_{1}n_{2})^{2}\log^{2}(n_{1}+n_{2})\right]$
$\displaystyle=\mathcal{O}\left[(n_{1}^{3}n_{2}^{2}+n_{1}^{2}n_{2}^{3})\log^{2}(n_{1}+n_{2})\right]$

∎

###  B Transcription Details

We split the handwritten text recognition task into three steps, those being
(1) page classification, (2) image segmentation, and (3) text recognition
(this is similar to the abstract approach described in [35]). For page
detection, we use an unsupervised approach based on clustering (specifically,
we use the DBSCAN algorithm of [47] on features extracted using a pre-trained
neural network). For image segmentation, we use a strategy based on point set
registration, the task of aligning points between an image and a template
[25]. We use the semantic segmentation model of [38] to extract vertical and
horizontal lines of the nurse records. We align the set of intersection points
of each page with a pre-specified template using the efficient probabilistic
point-set registration (FilterReg) by [54], thus obtaining a transformation
matrix for each image which we use, in combination with a pre-specified
overlay, to crop each field of interest into a separate image.

Having thus obtained segmented images of each field of interest, we use neural
networks to transcribe their contents. Here, we use a combination of
convolutional neural networks (CNN), those being based on the EfficientNetV2
architecture [98] and based on the approach in [37], and vision transformer
(ViT) architecture by [43]. We use different models depending on the field we
transcribe, as the content of the fields vary, with some consisting of
numbers, others characters, and yet others circling a number rather than
writing anything. Appendix Table B.1 provides information on all the different
groups of fields we consider. Note that some of these groups consist of many
fields with the same type of information (and as such the same “alphabet”,
consisting of the set of characters/digits that can occur in the field), while
others refer to just a single field. All fields are present on each nurse
record.

Table B.1: Grouping of Fields from the Nurse Records.

Group | #Fields | Maximum Length | Alphabet
---|---|---|---
breastfeed-7-do | 1 | 1 | $\\{1,2,3\\}^{\ast}$
dura-any-breastfeed | 1 | 2 | $\\{0,1,\dots,9\\}$
date | 7 | 4 | $\\{0,1,\dots,9\\}$
length | 2 | 3 | $\\{0,1,\dots,9\\}$
birth prior to due date | 1 | 1 | $\\{1,2\\}^{\ast}$
no. of weeks prior to due date | 1 | 2 | $\\{0,1,\dots,9\\}$
tab-b | 112 | 2 | $\\{0,1,\dots,9\\}$
weight | 8 | 5 | $\\{0,1,\dots,9\\}$
nurse-name (first) | 3 | $k^{\ast\ast}$ | $\\{a,b,\dots,\text{å}\\}$
nurse-name (last) | 3 | $k^{\ast\ast}$ | $\\{a,b,\dots,\text{å}\\}$

Notes: The table shows fields of the nurse records grouped together in such a
way that fields which are similar with respect to sequence length and alphabet
of their contents are put into one group. The first column refers to the name
given to the specific group of fields. The second column shows the number of
fields of the given group. The third column shows the maximum length of the
content of any field of the given group. The fourth column shows the alphabet
of the fields of the given group. ∗The alphabet of these fields are
specifically a circle being put around one of the digits shown, the digits
being pre-printed on the records. ∗∗While there is no clear limit to the
length of a name, the longest name in our training dataset contains 14
characters. Source: Table due to [12].

For each group of Appendix Table B.1, we train two neural networks, one based
on a CNN structure and one based on a ViT structure. Further, we create two
“meta”-groups, consisting of different groups for which a single alphabet
easily covers the different types of content. First, we create a
“Circle”-group consisting of breastfeeding at seven days old (breastfeed-7-do)
and born prior to due date (birth prior to due date), both of which consist of
circling a number. Second, we create an “Integer”-group consisting of duration
of breastfeeding (dura-any-breastfeed), length (length), number of weeks born
prior to due date (weeks prior to due date), Table B visit information
(tab-b), and weight (weight), as all groups consist of sequences of one to
five integers. For the “Cirlce”-group, we train both a CNN and a ViT, and for
the “Integer”-group we train only a ViT. In total, this results in training 23
neural networks, of which six models were chosen. For each group of Appendix
Table B.1, we evaluate the appropriate models and select the one that performs
best on a held-out test set not used for training. Appendix Table B.2 shows
information on the selected transcription models; for full details see [12]
and/or the official GitHub (https://github.com/TorbenSDJohansen/cihvr-
transcription, to be made public).

Table B.2: Model Differences – Selected Models.

Model |  Fields (Appendix Table B.1) | Resolution | Seq. len. | Batch size
---|---|---|---|---
| Panel A: ViT-based models
Circle |  breastfeed-7-do, birth prior to due date | 100x350 | 4 | 256
Integer seq. |  dura-any-breastfeed, length, weeks prior to due date, tab-b, weight | 90x230 | 7 | 1024
Last name |  nurse-name | 91x530 | 20 | 512
| Panel B: CNN-based models
Date |  date | 67x181 | 3 | 1024
First name |  nurse-name | 91x530 | 18 | 256
Weight |  weight | 80x258 | 5 | 1024

Notes: The table shows differences of hyperparameters of the final six neural
networks selected. The differences are all related to the groups of fields
used for training, as they differ in resolution (including aspect ratio) and
sequence length, which, together with varying availability of 1 vs. 2 GPUs,
led to different batch sizes. Note how the sequence lengths of this table
differ from those of Appendix Table B.1, often being longer. This is due to
the addition of certain special characters such as beginning of sequence and
end of sequence tokens being pre- and appended to the sequences, respectively,
for some models. Source: Table due to [12].

#### Transcription Performance

Appendix Table B.3 shows the transcription performance for the various fields
of the nurse records. These groups offer a slightly more detailed perspective
than those of Appendix Table B.1. Notably, the Table B group has been further
subdivided into 16 distinct groups, aligning with the 16 different columns
within the group, where each column corresponds to seven unique fields (visits
at months one, two, three, four, six, nine, and 12).

It is worth noting that we intentionally avoid character accuracy evaluation,
as our objective is to ensure the correctness of entire sequences. For
instance, in the context of character accuracy, transcribing a weight of 6,000
grams and transcribing 5,000 grams would yield a relatively high similarity
score, despite their significant practical difference. Our decision to
emphasize “sequence accuracy” helps mitigate such discrepancies.

Appendix Table B.3 also shows performance metrics that account for scenarios
where empty fields are excluded. This analysis demonstrates that our high
transcription accuracy is not solely based on accurately predicting a
substantial number of empty fields, which might be relatively straightforward
to achieve. Notably, the last column of the table indicates the proportion of
non-empty fields. It is important to note that this proportion does not
represent the share of non-empty fields across the entire collection of
records within a particular group. Instead, it specifically pertains to the
evaluation set used for testing purposes.

Table B.3: Automated Transcription Performance.

| Transcription accuracy (%) | Share non-empty (%)
---|---|---
| All | Non-empty |
Babbles | 92.9 | 97.5 | 61.5
Breastfeeding 7 days | 99.4 | 99.7 | 91.4
Care and cleanliness | 96.7 | 99.0 | 60.8
Date | 97.2 | 96.8 | 77.5
Duration breastfeeding | 97.6 | 97.9 | 97.3
Home economic status | 96.0 | 97.6 | 46.5
Home harmony | 97.7 | 98.9 | 26.5
In air | 91.7 | 97.5 | 61.2
Length | 99.0 | 99.0 | 97.2
Lifts head | 92.6 | 94.9 | 63.4
Mother daily hours working at home | 88.8 | 99.5 | 60.7
Mother daily hours working outside home | 93.0 | 99.5 | 70.8
Mother mental capacity | 96.8 | 97.4 | 40.9
Mother physical capacity | 97.1 | 98.1 | 41.6
Number of daily meals | 74.5 | 74.0 | 72.9
Nursery or kindergarten | 93.9 | 93.3 | 69.9
Nutrition | 96.2 | 97.0 | 80.2
Own bed | 96.0 | 99.7 | 59.7
Birth prior to due date | 99.0 | 99.5 | 88.4
Weeks prior to due date | 97.3 | 80.1∗ | 12.5
Sits | 91.9 | 91.8 | 70.3
Smiles | 93.0 | 98.0 | 61.5
Weight | 97.8 | 97.7 | 97.3
Nurse first name | 95.2 | 93.6 | 57.0
Nurse last name | 95.0 | 93.2 | 57.0

Notes: The table shows the accuracy (%) of the ML transcriptions for separate
groups of fields, measured on an independent test set not part of the data
used to train our neural networks. The second column shows the accuracy on the
full test sample. The third column shows the accuracy when excluding empty
fields. The fourth column shows the share of observations of the test set that
is non-empty for each group. ∗The low sequence accuracy for non-empty number
of weeks prior to due date is due to inconsistencies regarding manual
labelling of ranges such as “1-2”. In those cases, the label might either say
1 or 2, meaning that it is not possible to do better than guessing one of the
two numbers for a number of these cases. Allowing the number of weeks born
prior to the due date to differ by one increases the sequence accuracy to
95.3% for the non-empty cases and to 99.2% for the full sample. Source: Table
due to [12].

###  C Additional Figures and Tables

Figure C.1: Number of Children by Nurse District and Year of Birth.

Notes: The figure shows the number of children in our primary sample by nurse
district and year of birth.

Figure C.2: Rank-Order Correlation Between Mean Outcomes of Nurses.

Notes: The figure shows p-values and rank-order correlation coefficients from
Spearman rank-order correlation tests. For each variable, the mean value for
the children of each nurse is calculated, and the rank-order correlation
between nurses is then calculated (i.e., comparing the ranks of nurses for one
variable with the ranks of nurses for another variable). Given the number of
nurses (112), the p-values are calculated by means of a permutation test which
randomly drew 100,000 permutations and then calculated the share of
correlation coefficients below that of the non-permuted data, using this for
the p-value.

Figure C.3: Rank-Order Correlation Between Mean Outcomes of Nurse Districts.

Notes: The figure shows p-values and rank-order correlation coefficients from
Spearman rank-order correlation tests. For each variable, the mean value for
the children of each nurse district is calculated, and the rank-order
correlation between nurses districts is then calculated (i.e., comparing the
ranks of nurse districts for one variable with the ranks of nurse districts
for another variable). Given the number of nurse districts (16), the p-values
are calculated by means of a permutation test which randomly drew 100,000
permutations and then calculated the share of correlation coefficients below
that of the non-permuted data, using this for the p-value.

(a) 1959

(b) 1960

(c) 1961

(d) 1962

(e) 1963

(f) 1964

(g) 1965

(h) 1966

(i) 1967

Figure C.4: Box Plots of Average Years of Education by Nurses for Each Nurse
District and Year.

Notes: The figure shows box plots of average years of education of the
children by nurse for each district and year, meaning that each point
represents an average of the children allocated a specific nurse born in a
specific year. The different panels refer to different years (1959-1967) and
the different box plots within each panel refer to different nurse districts
(1-16).

(a) 1959

(b) 1960

(c) 1961

(d) 1962

(e) 1963

(f) 1964

(g) 1965

(h) 1966

(i) 1967

Figure C.5: Box Plots of Average Earnings During Ages 25-50 by Nurses for Each
Nurse District and Year.

Notes: The figure shows box plots of average earnings during ages 25-50 of the
children by nurse for each district and year, meaning that each point
represents an average of the children allocated a specific nurse born in a
specific year. The different panels refer to different years (1959-1967) and
the different box plots within each panel refer to different nurse districts
(1-16).

(a) 1959

(b) 1960

(c) 1961

(d) 1962

(e) 1963

(f) 1964

(g) 1965

(h) 1966

(i) 1967

Figure C.6: Average Years of Education (with Confidence Intervals) by Nurses
for Each Nurse District and Year.

Notes: The figure shows plots of average years of education with associated
95% confidence intervals (CIs) of the children by nurse for each district and
year, meaning that each point represents an average of the children allocated
a specific nurse born in a specific year. The different panels refer to
different years (1959-1967) and the different rows within each panel refer to
different nurse districts (1-16).

(a) 1959

(b) 1960

(c) 1961

(d) 1962

(e) 1963

(f) 1964

(g) 1965

(h) 1966

(i) 1967

Figure C.7: Average Earnings During Ages 25-50 (with Confidence Intervals) by
Nurses for Each Nurse District and Year.

Notes: The figure shows plots of average earnings during ages 25-50 with
associated 95% confidence intervals (CIs) of the children by nurse for each
district and year, meaning that each point represents an average of the
children allocated a specific nurse born in a specific year. The different
panels refer to different years (1959-1967) and the different rows within each
panel refer to different nurse districts (1-16).

(a) Breastfed at 1 mo.

(b) Breastfed at 6 mo.

(c) Duration of breastfeeding (mo.)

(d) Years of education

(e) Above mandatory education

(f) Share of time employed 25-50

(g) Average income 25-50

Figure C.8: Counterfactual Average Absolute Differences in Outcomes Between
Permuted Children of Different Nurses Within Same District-by-Year.

Notes: The figure shows histograms (with associated density plots) of averages
of absolute values of differences in average values of children with different
nurses but within same district-by-year group, but where children are permuted
in such a way as to be randomly allocated a nurse within the child’s district-
by-year group, with a total of 1000 permutations per panel. The black, solid
line in each panel shows the average of the underlying absolute differences in
the non-permuted sample, and the p-value reported in each panel is calculated
as the share of permutations that led to a more extreme value than the actual
average. The different panels refer to different outcome variables.

(a) Low BW

(b) Born prior to due date

(c) Weeks born prior to due date

(d) Born 1-3

(e) Female

(f) Father Missing

Figure C.9: Counterfactual Average Absolute Differences in Pre-Treatment
Variables Between Permuted Children of Different Nurses Within Same District-
by-Year.

Notes: The figure shows histograms (with associated density plots) of averages
of absolute values of differences in average values of children with different
nurses but within same district-by-year group, but where children are permuted
in such a way as to be randomly allocated a nurse within the child’s district-
by-year group, with a total of 1000 permutations per panel. The black, solid
line in each panel shows the average of the underlying absolute differences in
the non-permuted sample, and the p-value reported in each panel is calculated
as the share of permutations that led to a more extreme value than the actual
average. The different panels refer to different pre-treatment variables.


# Data-driven Invariance for Reference Governors

Ali Kashani<EMAIL_ADDRESS>Claus Danielson<EMAIL_ADDRESS>Department of
Mechanical Engineering, University of New Mexico, Albuquerque, NM, US 87106

###### Abstract

This paper presents a novel approach to synthesizing positive invariant sets
for unmodeled nonlinear systems using direct data-driven techniques. The data-
driven invariant sets are used to design a data-driven reference governor that
selects a reference for the closed-loop system to enforce constraints. Using
kernel-basis functions, we solve a semi-definite program to learn a sum-of-
squares Lyapunov-like function whose unity level-set is a constraint
admissible positive invariant set, which determines the constraint admissible
states as well as reference inputs. Leveraging Lipschitz properties of the
system, we prove that tightening the model-based design ensures robustness of
the data-driven invariant set to the inherent plant uncertainty in a data-
driven framework. To mitigate the curse-of-dimensionality, we repose the semi-
definite program into a linear program. We validate our approach through two
examples: First, we present an illustrative example where we can analytically
compute the maximum positive invariant set and compare with the presented
data-driven invariant set. Second, we present a practical autonomous driving
scenario to demonstrate the utility of the presented method for nonlinear
systems.

††thanks: This material is based upon work supported by the National Science
Foundation under NSF Grants CMMI-2105631 and CMMI-2303157. Any opinions,
findings, and conclusions or recommendations expressed in this material are
those of the authors and do not necessarily reflect the views of the National
Science Foundation.

,

##  1 Introduction

Safety-critical cyber-physical systems, such as autonomous vehicles, surgical
systems, and industrial robots, must operate within strict constraints to
avoid damaging the machine or its environment, including the human operators.
These constraints are highly application-oriented. Thus, the existing
controller may not be primarily designed for such constraints and be only
responsible for general issues like stability and reference tracking. For
example, the primary controller in an assisted-driving vehicle may be designed
for stability, while driver’s maneuvers can cause the car to cross lanes. One
solution is creating a constraint-enforcing controller (e.g. model predictive
control (mpc)) to replace the existing controller. This is often undesirable
since it discards the careful testing and tuning of the existing controller.
Furthermore, control design techniques such as mpc typically require a model
which may not be available. In the assisted driving example, an uncertain
model may lead to impulsive maneuvers that destabilize the system e.g., cause
the car to roll, or at least discomfort the passengers [1].

Alternatively, a reference governor (rg) [2] can be added to the existing
closed-loop system as an outer-loop constraint-enforcing controller that
adjusts the reference commanded to the existing controller. This is
advantageous since redesigning the existing controller for every new scenario
is challenging, while rg is implemented as an add-on to accommodate new
scenarios. Moreover, the proposed direct data-driven rgs can be designed
without having the information about the system that was required for
designing the existing controller.

Solving rgs can be intractable when the system is mathematically inaccessible,
complex and uncertain. Some formulations of rgs (e.g. [3]) require infinite
horizon planning to ensure the current reference is safe, which can be
mitigated by using high-level dynamical properties such as invariance [4]. In
this context, a controlled invariant (ci) set is a subset of the state-space
for which there exists a reference trajectory such that the state trajectory
does not leave that set. The ci sets reside on the existence of a reference
that ensures invariance. However, computing the feasible references is not
trivial. Even if the plant model and a ci set are mathematically accessible,
computing the reference may be intractable e.g., require solving a non-convex
optimization problem. This is more challenging in a data-driven approach. For
data-sets of the form $\\{x_{t+1},x_{t},r_{t}\\}$, we only have one reference
$r_{t}$ associate with each state $x_{t}$ at each time $t$, which may or may
not promote invariance since the ci set is not known in advance.

We simplify the problem of computing a ci set to the computation of a
collection of pi sets for different references whose union is a ci set. Once
the pi sets are computed, they form an admissible set that reduces the rg to a
one step planner; find a reference that keeps the state inside the admissible
set. When the model of the system is available, the invariant sets can be
computed by standard algorithms as in [5]. These algorithms can compute
maximal output admissible sets for linear systems as in [6], that promotes
constraint enforcement. This form of rg is popular for linear systems (e.g.
[7, 8, 9]) as well as nonlinear systems (e.g., [10]). However, when the model
is not available, these algorithms require indirect data-driven approach
wherein a model is first identified from data and then used for model-based
controller synthesis. Moreover, computing the maximal invariant sets may be
overly ambitious. To the best of our knowledge, a general data-driven approach
that finds the exact maximal invariant sets does not exist. Even the model-
based standard algorithms [5], can be computationally intractable. Thus, we
will settle for a sub-maximal invariant set, provided we can guarantee its
invariance despite the uncertainty due to data-driven framework. This is
challenging since approximations of an invariant set are not approximately
invariant, e.g. [11]. Thus, an approximate invariant set may allow constraint
violates to occur.

This paper present an alternative direct data-driven approach wherein the
invariant set is learned directly from data. This approach is attractive since
it skips the system identification step. A popular method is parameterizing an
invariant set as the level-set of a characterizing function. The data-driven
computation of the characterizing function makes it possible to leverage
machine learning tools. Learning-based techniques show promise in constructing
Lyapunov-like and barrier certificates from data that support invariance.
Neural networks that learn Lyapunov function of the system, can be used to
compute either invariance for autonomous systems (e.g. [12, 13]) or design
stabilizing control policies (e.g. [14]). The learning process usually
includes solving a semi-definite program (sdp), which is sometimes dealt with
sum of squares (sos) programming for polynomial systems with polynomial
constraints [15, 16, 17]. Our approach is inspired by sos programming but is
not limited to polynomial systems. A similar approach is introduced in [18]
though they do not discuss how to adjust the control or reference policy, but
to compute positive invariance for a certain policy. Another similar approach
in [19] includes the control design but requires a linear model that
represents the nonlinear system at steady state. There also data-driven rg
methods (e.g. [3]) which do not use invariance principle and try to directly
solve a finite-time receding horizon optimization problem using machine
learning tools, for which computational cost and tractability as well as
formal safety guarantees need to be investigated.

This paper employs learning-based techniques to develop data-driven reference
governors for nonlinear systems. We outline the necessary information about
the plant in section 2. Then, we illustrate the transformation of the problem
of identifying invariant sets into the creation of a Lyapunov-like function,
and detail how to learn this function using either the model of the system or
available data in section 3. Proofs are provided to ensure robustness against
uncertainties due to data-driven frameworks and equilibrium approximations. In
section 4, we present two illustrative examples: firstly, a standard second-
order linear system to showcase our method, and secondly, the nonlinear lane-
keeping problem to demonstrate the effectiveness of our approach.

#### Notations and Definitions

The Frobenius inner product of square matrices $A,B$ is denoted by $\langle
A,B\rangle=\text{Tr}(A^{\top\\!}B)$, where Tr$(.)$ is the trace operator.
Matrix inequality $A\succ B$ means $A-B\succ 0$ is positive definite. The
successor of the state $x_{t}$ is denoted by $x^{+}=x_{t+1}$. The argument(s)
of the functions and time index $t$ are sometimes dropped to simplify the
notation. The Cartesian product and the relative complement of a set are
denoted by $\times$ and $\backslash$ operators, respectively.

For an autonomous system $x^{+}=f(x)$, a set $\mathcal{O}$ is pi if for any
$x\in\mathcal{O}$ we have $f(x)\in\mathcal{O}$. For a controlled system
$x^{+}=f(x,r)$, a set $\mathcal{C}$ is ci if for any $x\in\mathcal{C}$ there
exists a reference $r\in\mathcal{R}$ such that $f(x,r)\in\mathcal{C}$. The
maximal pi set $\mathcal{O}_{\infty}$ (maximal ci set $\mathcal{C}_{\infty}$)
is the pi (ci) set which contains all other pi (ci) sets
$\mathcal{O}\subseteq\mathcal{O}_{\infty}$
($\mathcal{C}\subseteq\mathcal{C}_{\infty}$) that are subsets of the
constraint set $\mathcal{O}\subseteq\mathcal{X}$
($\mathcal{C}\subseteq\mathcal{X}$).

##  2 Reference Governor Synthesis Problem

Our objective is to use data to achieve constrained reference tracking for an
unmodeled system. Our knowledge of this system is limited to a data-set
comprised of state trajectories corresponding to a reference $r\in\mathbb{R}$.
Specifically, every reference
$r=\bar{r}\in\\{\bar{r}^{i}\\}_{i=1}^{N_{\bar{r}}}$ has generated $N_{T}$
state trajectories $\textbf{x}_{j}^{i}\in\mathbb{R}^{n}\times[0,T]$ of length
$T$ in the data-set

$\mathcal{D}=\\{\\{\textbf{x}_{j}^{i}\\}_{j=1}^{N_{T}},~{}\bar{r}^{i}\\}_{i=1}^{N_{\bar{r}}},$
(1)

where the length $T$ and number $N_{\bar{r}}$ of the trajectories are
considered constant for notational simplicity. For the purpose of proving
invariance, we assume that the data (1) comes from an unmodeled discrete-time
dynamical system

$x^{+}=f(x,r),\quad x\in\mathcal{X},~{}r\in\mathcal{R},\\\ $ (2)

where $r\in\mathcal{R}$ is the reference in the feasible range
$\mathcal{R}\subseteq\mathbb{R}$, and $x\in\mathcal{X}\subseteq\mathbb{R}^{n}$
is the state that is aimed to stay in the state constraint set

$\mathcal{X}:=\\{x|~{}g(x)\leq 1\\},$ (3)

which is expressed as the sub-level set of some Lipschitz continuous function
$g$. We pose the following assumptions about the unmodeled system (2).

######  Assumption 1.

(System Assumptions)

  1. (a)

The state $x$ is available, but the dynamics $f$ are unmodeled.

  2. (b)

The system (2) is asymptotically input-to-state stable (iss).

  3. (c)

The dynamics $f$ are Lipschitz continuous with known Lipshitz constant
$L_{f}$.

Assumption 1(a) generally holds for dynamical systems where the state $x$ is
measurable. Alternatively, we can define the state $x$ to be a window of any
_flat_ measurement $\\{y_{k}\\}_{k=t-n+1}^{t}$ about the system when the
unmodeled dynamics $f$ is _differentially flat_ [20]. Assumption 1(b) ensures
that for constant reference input $r_{t}=\bar{r}$, the state converges to a
unique equilibrium $x_{t}\to x_{\infty}(\bar{r})$ as $t\to\infty$. This holds,
for example, when the system is equipped with an inner-controller. There are
several reasons that the closed-loop dynamics (2) is considered unmodeled
despite the existence of the inner-controller:

  * •

The inner-controller was designed using model-free techniques e.g. a hand-
tuned PID controller.

  * •

The plant was manufactured by a third party that is unwilling to share their
proprietary models with us.

  * •

The inner-controller is a complicated collection of legacy controllers and
switching logic.

  * •

A high-fidelity model is available but is too high-dimensional or complex for
analytic computation.

Assumption 1(c) will be used to guarantee robustness to plant uncertainty due
to the data-driven framework, since we have finite-data (1) rather than a
high-fidelity model (2). We note that all items of Assumption 1 are consistent
with the literature [19].

###  2.1 Reference Governors

rgInner-ControllerPlant$r_{d}$$r$$u$$x$Inner-loopOuter-loop$x^{+}=f(x,r)$
Figure 1: Reference governor (rg) design for a system equipped with an inner-
controller.

An rg is an outer-loop controller, as shown in Fig. 1, that minimally adjusts
the reference $r$ from the desired reference $r_{d}$ to ensure the constraints
$\mathcal{X}$ of the system (2) are not violated $x\in\mathcal{X}$. Often this
is posed as a dynamic constrained optimization problem with a sufficiently
long prediction horizon where the dynamics $f$ are known in order to predict.
Alternatively, an rg can be reduced to a static optimization problem, if an
admissible set of the system (2) is available. Then, rg can be posed as

$\displaystyle r^{\star}\in\arg\min\>$ $\displaystyle\|r-r_{d}\|$ (4a) s.t.
$\displaystyle\begin{bmatrix}x\\\ r\end{bmatrix}\in\Lambda,$ (4b) where
$\Lambda$ is an admissible set
$\Lambda\\!=\\!\\{(x,r)|~{}r\\!\in\\!\mathcal{R},~{}f(x,r)\in\mathcal{C}\\},$
(4c)

which is derived from a ci set $\mathcal{C}\subseteq\mathcal{X}$ of the system
(2).

Constraint enforcement is guaranteed when rg (4) is recursively feasible
$x_{t}\in\mathcal{C}\subseteq\mathcal{X},\forall t$. The following lemma shows
that the existence of a ci set is essential for recursive feasibility of rg
(4).

######  Lemma 1.

rg (4) is recursively feasible for all the states
$x\in\textnormal{Pre}(\mathcal{C})$ in the one-step predecessor
$\textnormal{Pre}(\mathcal{C})$ of $\mathcal{C}$ if and only if $\mathcal{C}$
is ci.

Since $x\in\textnormal{Pre}(\mathcal{C})$, by definition there exist a
reference input that takes the state successor $f(x,r)$ to
$\mathcal{C}\subseteq\mathcal{X}$. If $\mathcal{C}$ is ci, then there exist a
reference input at any time that keeps the state in
$\mathcal{C}\subseteq\mathcal{X}$. If $\mathcal{C}$ is not ci, then there
exists a point $x\in\mathcal{C}$ that state trajectory starting from $x$ will
eventually leave $\mathcal{C}$. Lemma 1 suggests that the sufficient plant
knowledge required for synthesizing a recursively feasible rg (4) is only an
admissible set $\Lambda$ (4c), which includes a ci set
$\mathcal{C}\in\mathcal{X}$ and the range of the reference $r$ for every point
in $\mathcal{C}$ that renders $\mathcal{C}$ being ci.

The maximal ci set $\mathcal{C}_{\infty}\subseteq\mathcal{X}$ is the largest
set in the constraint set $\mathcal{X}$ that renders rg recursively feasible
by Lemma 1 since
Pre$(\mathcal{C}_{\infty})\cap\mathcal{X}=\mathcal{C}_{\infty}$. Consequently,
rg (4) has the largest domain when the maximal ci set $\mathcal{C}_{\infty}$
is used, which produces the maximal admissible set $\Lambda_{\infty}$ of the
system (2) for (4c).

In the next section we will show that our approximation
$\Lambda\subseteq\Lambda_{\infty}$ is an admissible set since it will be
derived from a sub-maximal ci set $\mathcal{C}\subseteq\mathcal{C}_{\infty}$.
We note that finding the feasible reference $r$ that satisfy
$f(x,r)\in\mathcal{C}$ may not be tractable even if the model $f$ and the ci
set $\mathcal{C}$ are known e.g., (4) may be a non-convex optimization. In
section 3, we will present an algorithm that provides a tractable
implementation to rg (4).

##  3 Data-driven Reference Governor

In this section, first we describe a model-based method to compute a sub-
maximal ci set. Then, we extend the results to a data-driven approach.

###  3.1 PI Inner-Approximations of Maximal CI

The definition of a ci set requires the existance of a reference $r$ such that
$f(x,r)\in\mathcal{C}$. This presents an issue for data-driven invariance
since our data-set may not contain the appropriate references. In this
section, we overcome this issue by inner-approximating a ci set using a union
of pi sets which are more amenable to data-driven synthesis.

We use the fact that once the reference policy $r_{t}=\bar{r}$ is determined
for all time $t$, the system (2) becomes an autonomous system with a
particular dynamics $f(x;\bar{r})$ parameterized by the reference policy
$\bar{r}$. Now we can define pi sets for the autonomous dynamics
$f(x;\bar{r})$, and use the pi sets to find a ci set for the system (2). Every
pi set $\mathcal{O}({\bar{r}})\in\mathcal{X}$ of the autonomous dynamics
$f(x;\bar{r})$ parameterized by the reference policy $\bar{r}\in\mathcal{R}$
is a subset of the maximal ci set $\mathcal{C}_{\infty}$ of the system (2).
The following proposition shows that the union

$\mathcal{C}=\bigcup_{\bar{r}\in\bar{\mathcal{R}}}\mathcal{O}({\bar{r}})\subseteq\mathcal{C}_{\infty}\subseteq\mathcal{X}$
(5)

of the resulting pi sets $\mathcal{O}({\bar{r}})$ over any subset
$\bar{\mathcal{R}}\subseteq\mathcal{R}$ is a ci set.

######  Proposition 1.

Let $\mathcal{O}(\bar{r})$ be pi for the autonomous dynamics $f(x;\bar{r})$
where $\bar{r}\in\bar{\mathcal{R}}$. Then, (5) is ci for the system (2).

At the state $x\in\mathcal{C}$ we have $x\in\mathcal{O}(\bar{r})$ for some
$\bar{r}\in\bar{\mathcal{R}}$. Therefore, the reference input $\bar{r}$ keeps
the state successor $x^{+}\in\mathcal{O}(\bar{r})$, thus renders (5) ci. The
advantage of (5) is that it simplifies the problem of computing a ci set
$\mathcal{C}$ to the computation of pi sets $\mathcal{O}({\bar{r}})$. This is
especially helpful in data-driven approach since ci sets consider the
existance of a reference that ensures invariance. Conversely, the ci set (5)
includes all admissible references supported by data, where the corresponding
admissible state regions are described by the pi sets
$\mathcal{O}({\bar{r}})$.

As a result of Proposition 1, rg (4) can be revised by

$\displaystyle r^{\star}\in\arg\min\>$ $\displaystyle\|\bar{r}-r_{d}\|$ (6a)
s.t. $\displaystyle\begin{bmatrix}x\\\ \bar{r}\end{bmatrix}\in\bar{\Lambda},$
(6b) where the admissible set
$\bar{\Lambda}=\\{(x,\bar{r}):\bar{r}\in\bar{\mathcal{R}}~{}\text{s.t.}~{}x\in\mathcal{O}({\bar{r}})\\}$
(6c)

enables the synthesis of rg by having the pi sets $\mathcal{O}({\bar{r}})$ for
every reference policy ${\bar{r}}\in\bar{\mathcal{R}}$. The new representation
(6) of rg is now an integer program (ip), which selects between the policies
${\bar{r}}$ that enforces constraints $\mathcal{X}$. Algorithm 1 provides a
tractable solution to this ip (6), by searching for the nearest reference
policy $\bar{r}\in\bar{\mathcal{R}}$ such that $x_{t}\in\mathcal{O}(\bar{r})$
the current state $x_{t}$ lies in the corresponding pi set
$\mathcal{O}(\bar{r})$. Moreover, recursive feasibility of rg (6) is verified
in the following corollary.

######  Corollary 1.

The rg (6) is recursively feasible if $\mathcal{O}({\bar{r}})$ is pi for every
reference policy ${\bar{r}}\in\bar{\mathcal{R}}$.

The proof follows the proof of Lemma 1, where
$\mathcal{C}\subseteq\mathcal{X}$ is ci by Proposition 1, and
$x\in\mathcal{C}$ is guaranteed for all time since $\mathcal{O}(\bar{r})$ is
pi.

Algorithm 1 Solving rg (6)

1:input: current state $x$, desired reference $r_{d}$, and
$\\{\mathcal{O}(\bar{r}^{i})\\}_{i=1}^{N_{\bar{r}}}$ (or
$\\{\mathcal{D}_{\mathcal{O}}^{i}\\}_{i=1}^{N_{\bar{r}}}$ from Alg. 2)

2:initiate:
$i^{\star}=\arg\min_{\bar{r}^{i}\in\bar{\mathcal{R}}}\>\|\bar{r}^{i}-r_{d}\|$

3:while $x\not\in\mathcal{O}({\bar{r}^{i^{\star}}})$ do

4: $\bar{\mathcal{R}}=\bar{\mathcal{R}}\backslash\\{\bar{r}^{i^{\star}}\\}$

5:
${i^{\star}}=\arg\min_{\bar{r}^{i}\in\bar{\mathcal{R}}}\>\|\bar{r}^{i}-r_{d}\|$

6:end while

7:output: $\bar{r}^{i^{\star}}$

###  3.2 Model-Driven Approach

In this section, we present a model-based optimization algorithm for computing
the pi sets $\mathcal{O}(\bar{r})$. The sequel will adapt this approach for a
data-driven paradigm.

We use sub-level representation
$\mathcal{O}({\bar{r}})\\!=\\!\\{x|~{}V(x;\bar{r})\\!\leq\\!1\\}$ to transform
the problem of synthesizing a set $\mathcal{O}({\bar{r}})\subseteq\mathcal{X}$
into the problem of synthesizing a function $V(x;\bar{r})$. This allows us to
leverage machine learning techniques such as neural barrier/Lyapunov functions
[21, 14] to estimate a pi set $\mathcal{O}(\bar{r})$. The conditions that the
function $V$ must satisfy are

$\displaystyle x\not\in\mathcal{X}$ $\displaystyle\mapsto V(x~{};\bar{r})\geq
1,$ (7a) $\displaystyle x\in\mathcal{O}(\bar{r}):V(x~{};\bar{r})\leq 1$
$\displaystyle\mapsto V(f(x;\bar{r})~{};\bar{r})\leq 1,$ (7b)

where (7a) ensures $\mathcal{O}\subseteq\mathcal{X}$ by the S-procedure, and
(7b) ensures the trajectories starting in $\mathcal{O}$ do not leave
$\mathcal{O}$. However, satisfying (7) is not tractable since (7b) must hold
inside the unknown set $\mathcal{O}({\bar{r}})$. Thus, similar to neural
barrier/Lyapunov functions approaches, we relax the condition (7) with

$\displaystyle x\not\in\mathcal{X}$ $\displaystyle\mapsto V(x~{};\bar{r})\geq
1,$ (8a) $\displaystyle x\in\mathcal{X}$ $\displaystyle\mapsto
V(f(x;\bar{r})~{};\bar{r})\leq(1\\!-\\!\gamma)~{}V\\!(x~{};\bar{r})+\gamma\\!-\\!\beta,$
(8b)

where the Lyapunov-like condition (8b) implies (7b). The parameter
$0\leq\gamma\leq 1$ forms a convex combination of $V$ and $1$ in the right
hand side of (8b) that dictates $V$ to be decreasing outside $\mathcal{O}$.
However, it allows $V$ to be increasing inside $\mathcal{O}$, similar to a
barrier function. The drift parameter $0\leq\beta\leq\gamma$ enforces $V$ to
be decreasing in the region
$\\{x|~{}V\\!\geq\\!\frac{\gamma-\beta}{\gamma}\\}$ that includes a strip
inside and close to the boundary of $\mathcal{O}$. Although invariance is
known to be fragile at the boundary, a nonzero $\beta$ can provide robustness
to the uncertainty of data-driven framework.

We define the pi set $\mathcal{O}(\bar{r})\subseteq\mathcal{X}$ of the
autonomous dynamics $f(x;\bar{r}$), as a sub-level set of the Lyapunov-like
function

$V(x~{};\bar{r})=\varphi(x;\bar{r})^{\top\\!}P\varphi(x;\bar{r})=\langle
P,\varphi(x;\bar{r})\varphi(x;\bar{r})^{\top\\!}\rangle,$ (9)

where the positive definite matrix $P\in\mathbb{R}^{N_{\phi}\times
N_{\phi}}\succ 0$ parameterizes the set $\mathcal{O}(\bar{r})$ and,

$\varphi(x~{};\bar{r})=\phi(x)-\phi(x_{\infty}(\bar{r})),$ (10)

and $\phi(x)\in\mathbb{R}^{N}_{\phi}$ consists of $N_{\phi}$ Lipschitz
continuous scalar dictionary functions

$\phi(x)=[\phi_{1}(x),\cdots,\phi_{N_{\phi}}(x)]^{\top\\!}.$

For the pi sets $\mathcal{O}(\bar{r})$ in rg (6c), we consider the following

$\mathcal{O}({\bar{r}})\\!:=\\!\\{x|~{}V\\!\leq\\!\big{(}1-c^{\top\\!}\phi(x_{\infty}(\bar{r}))\big{)}^{2}\\!\\!,~{}\text{and}~{}c^{\top\\!}\phi(x_{\infty}(\bar{r}))\\!\leq\\!1\\},$
(11)

where $c\in\mathbb{R}^{N_{\phi}}$ parameterizes an inner-approximation
$\hat{\mathcal{X}}:=\\{x|~{}c^{\top\\!}\phi(x)\leq 1\\}\subseteq\mathcal{X}$
of the constraint set $\mathcal{X}$. The argument
$c^{\top\\!}\phi(x_{\infty}(\bar{r}))\\!\leq\\!1$ simply means the equilibrium
is admissible $x_{\infty}\in\hat{\mathcal{X}}$. To obtain $c$, one can solve
the following linear program (lp)

$\displaystyle c=\arg\min\>$
$\displaystyle\int_{\mathcal{X}}c^{\top\\!}\phi~{}\mathrm{d}\delta,$ (12a)
s.t. $\displaystyle c^{\top\\!}\phi(x)\geq g(x),~{}~{}\forall x$ (12b)

where $g(x)$ characterizes the constraint set $\mathcal{X}$ from (3), and the
cost (12a) is a Lebesgue integral with some measure $\delta$ over
$\mathcal{X}$. However, we can avoid solving the lp (12) by adding $g(x)$ as
an element of $\phi(x)$. For instance, $c=[1,0,\cdots,0]^{\top\\!}$ is the
trivial solution to (12) when $\phi_{1}(x)=g(x)$. We note that this does not
necessarily increase the volume of the pi set
$\mathcal{O}(\bar{r})\subseteq\hat{\mathcal{X}}\subseteq\mathcal{X}$, since
the choice of $\phi$ effects the Lyapunov-like function.

Before presenting our data-driven approach, let us consider the case where the
autonomous dynamics $f(x;\bar{r})$, and thus, the corresponding equilibrium
$x_{\infty}(\bar{r})$ are known. In this case, the following sos program may
be solved to find the parameterizing positive definite matrix $P$,

$\displaystyle P\\!\in\\!\arg\min$ $\displaystyle\langle
P^{\prime},\Psi\rangle$ (13a) s.t. $\displaystyle\langle
P^{\prime},(\varphi\\!\circ\\!f)(\varphi\\!\circ\\!f)\\!^{\top\\!}\\!\\!-\\!(1\\!-\\!\gamma)\varphi\varphi\\!^{\top\\!}\rangle\\!\leq\\!\gamma,~{}\forall
x$ (13b) $\displaystyle cc^{\top\\!}\preceq P^{\prime},$ (13c)

where $\varphi\circ f=\varphi(f(x;\bar{r})~{};\bar{r})$, and

$\Psi=\int_{\mathcal{X}}\\!\varphi\varphi^{\top\\!}~{}\mathrm{d}\delta.$ (14)

The sos program (13) minimizes a Lebesgue integral of the Lyapunov-like
function (9) with some measure $\delta$ over $\mathcal{X}$, since ideally we
like to find maximal pi sets for (11) to derive $\mathcal{C}$ (5) that
produces a sub-maximal ci set. However, the choice of the lift $\phi$ plays
the key role in enforcing (13b). Classic sos programming uses polynomials
basis functions $\phi$ where the dynamics $f$ belong to the space of
polynomials [17]. Kernel methods in sos programming use functional polynomials
to capture nonlinearities, where the dynamics become polynomials after a
nonlinear transformation [22]. In our numerical results, we use machine-
learning kernel functions to train a sos Lyapunov-like neural network [18].

###  3.3 Data-Driven Approach

In this section, we adapt the model-based pi set synthesis problem (13) to a
data-driven paradigm.

Algorithm 2 summarizes our data-driven method for calculating the sub-maximal
ci set (5). It uses the data set (1) to compute the sub-maximal pi sets (11)
for different reference policies $\bar{r}\in\bar{\mathcal{R}}$. This
information can be used to solve rg (6) by Algorithm 1, when the dynamics $f$
are not modeled. The main operation in Algorithm 2 is computing $P$ from the
data (1). This requires transforming the model-based method (13) into a data-
driven method, which is presents two challenges.

First, the equilibrium point $x_{\infty}(\bar{r})$ is not necessarily
available. Thus, we use an ensemble average to approximate

$x_{\infty}(\bar{r}^{i})\approx\frac{1}{N_{T}}\sum_{j=1}^{N_{T}}{\textbf{x}_{j}^{i}}_{T}$
(15)

for every $\bar{r}\in\\{\bar{r}^{i}\\}_{i=1}^{N_{\bar{r}}}$, where
${\textbf{x}_{j}^{i}}_{T}$ denotes the last point on the $j$th trajectory from
the data set $\eqref{eq:data}$ that is generated by the reference policy
$r_{t}=\bar{r}^{i}$. For the approximation (15) to be reliable, every
reference policy $\bar{r}$ should asymptotically take the system to the unique
equilibrium $x_{T}\to x_{\infty}(\bar{r})$ as $T\to\infty$. This is achievable
by Assumption 1(b), since we can generate converging trajectories using
constant values of $\bar{r}$. Later in Theorem 2, we will prove that our pi
sets are robust to the approximation (15) with finite $T$.

Second, since we only have finite data, our plant is uncertain; The condition
(13b) must hold for every point $x\in\mathcal{X}$, while in a data-driven
paradigm only a finite number $N_{s}$ of sample pairs
$\\{(x_{k}^{+},x_{k})\\}_{k=1}^{N_{s}}$ are available. Thus, to insure the pi
set $\mathcal{O}(\bar{r})$ (11) is robust to plant uncertainty due to the
data-driven framework i.e., enforcing (13b) on finite data follows (13b) holds
everywhere in the pi set $\mathcal{O}(\bar{r})$, we tighten (13b) by the term

$\epsilon_{k}\geq
2L_{\phi}\delta(L_{f}\|\varphi_{k}^{+}\|+\|\varphi_{k}\|)+(L_{\phi}L_{f}\delta)^{2}$
(16)

to the Lyapunov-like constraint (13b), where $L_{\phi}$ and $L_{f}$ are the
Lipschitz constants of the dictionary functions $\phi$ and dynamics $f$,
respectively, and $\delta$ is the density of the data samples i.e., for any
$x\in\mathcal{X}$ there exists a sample $x_{k}$ such that
$\|x-x_{k}\|\leq\delta$. We note that the $\delta$ and $L_{\phi}$ are
available; $\delta$ can be computed from the data-set and $L_{\phi}$ depends
on our choice of dictionary functions. The Lipschitz constant $L_{f}$ can be
computed from data using e.g. [23]. The function $\varphi$ from (10) is used
to “lift” the sample pairs $\\{(x_{k}^{+},x_{k})\\}_{k=1}^{N_{s}}$,

$\varphi_{k}=\varphi(x_{k};\bar{r}),~{}~{}~{}\varphi_{k}^{+}=\varphi(x_{k}^{+};\bar{r}).$
(17)

Our data-driven approach finds the matrix $P$ for the sub-maximal pi set
$\mathcal{O}(\bar{r})$ (11) by the following sdp

$\displaystyle P\in\arg\min~{}$ $\displaystyle\langle P^{\prime},\Psi\rangle$
(18a) s.t. $\displaystyle\langle
P^{\prime},\psi_{k}\rangle\leq\gamma-\epsilon_{k}\lambda,~{}k\\!=\\!1,\\!\cdots\\!,\\!N_{s},$
(18b) $\displaystyle cc^{\top\\!}\preceq P^{\prime}\preceq\lambda
I_{N_{\phi}},$ (18c)

where
$\psi_{k}=\varphi_{k}^{+}\varphi_{k}^{+{\top\\!}}-(1-\gamma)\varphi_{k}\varphi_{k}^{\top\\!}$,
and $\epsilon_{k}$ satisfies (16). We may either treat $\lambda$ as a decision
variable or pre-assign it in sdp (18). Furthermore, $\Psi$ in the cost (18a)
can be analytically calculated by (14) since $\varphi$ is available in the
analytical form (10). Alternatively, we can replace the integral $\Psi$ (14)
with a summation over the outer product of the lifted sample pairs (17),

$\Psi=\sum_{k=1}^{N_{s}}\varphi_{k}\varphi_{k}^{\top\\!}\delta(x_{k}),$

with some measure $\delta$.

Algorithm 2 Data-driven sub-maximal ci set computation

1:input: Data set $\mathcal{D}$ from (1).

2:for $i=1,\cdots,N_{\bar{r}}$ do

3: approximate the equilibrium $x_{\infty}(\bar{r}^{i})$ by (15)

4: for $k=1,\cdots,N_{s}$, where $x_{k}\in\mathcal{D}^{i}$ do

5: lift  $\displaystyle\varphi_{k}$
$\displaystyle=\phi(x_{k})-\phi(x_{\infty}(\bar{r}^{i}))$
$\displaystyle\varphi_{k}^{+}$
$\displaystyle=\phi(x_{k}^{+})-\phi(x_{\infty}(\bar{r}^{i}))$

6: end for

7: solve sdp (18) (or lp (22) then (21)) for $P$,

8: store
$\mathcal{D}_{\mathcal{O}}^{i}:=\\{x_{\infty}(\bar{r}^{i}),P^{i}=P\\}$,

9:end for

10:output: $\\{\mathcal{D}_{\mathcal{O}}^{i}\\}_{i=1}^{N_{\bar{r}}}$

The following theorem shows that the solution $P$ to sdp (18), if feasible,
produces a pi set $\mathcal{O}({\bar{r}})\subseteq\mathcal{X}$ (11) for the
autonomous dynamics $f(x;\bar{r})$. We will later demonstrate the feasibility
of sdp (18).

######  Theorem 1.

Let Assumption 1 hold. Let $c$ and $P$ be the solutions to (12) and (18),
respectively. Let the data samples have a density of $\|x-x_{k}\|\leq\delta$.
Then, $\mathcal{O}({\bar{r}})$ from (11) is constraint admissible
$\mathcal{O}({\bar{r}})\subseteq\mathcal{X}$, and pi for the autonomous
dynamics $f(x;\bar{r})$ of the system (2) with reference policy
$r_{t}=\bar{r}$.

First, we show $\mathcal{O}(\bar{r})$ (11) is constraint admissible
$\mathcal{O}({\bar{r}})\in\mathcal{X}$. The lower bound (18c) on $P$ enforces
the scalar inequality
$V\geq(c^{\top\\!}(\phi-\phi(x_{\infty}(\bar{r}))))^{2}$, which by the
S-procedure it follows

$\mathcal{O}({\bar{r}})\subseteq\\{x|~{}c^{\top\\!}(\phi-\phi(x_{\infty}(\bar{r})))\leq(1-c^{\top\\!}\phi(x_{\infty}(\bar{r})))\\}.$

Again by the S-procedure, we get

$\mathcal{O}({\bar{r}})\subseteq\\{x|~{}c^{\top\\!}\phi\leq
1\\}\subseteq\mathcal{X}$

since $c^{\top\\!}\phi(x)\geq g(x)$ is enforced by (12b).

Next, we will show $\mathcal{O}({\bar{r}})$ is pi. The Lyapunov-like function
(9) is purely quadratic $V=\langle
P,\varphi\varphi^{\top\\!}\rangle=\|\varphi\|_{P}^{2}$ with respect to
$\varphi$, where $\|\varphi\|_{P}=\varphi^{\top\\!}P\varphi$ is the weighted
2-norm. Therefore, the Lyapunov-like function $V$ (9) can be upper-bounded at
any state successor $x^{+}$ with respect to an adjacent successor sample
$x_{k}^{+}$,

$\|\varphi^{+}\|_{P}^{2}\leq\|\varphi_{k}^{+}\|_{P}^{2}+2\|\varphi^{+}-\varphi_{k}^{+}\|_{P}\|\varphi_{k}^{+}\|_{P}+\|\varphi^{+}-\varphi_{k}^{+}\|_{P}^{2},$

where $\varphi^{+}=\varphi(x^{+};\bar{r})$. By the Lipschitz continuity, we
can conclude

$\|\varphi^{+}\|_{P}^{2}\leq\|\varphi_{k}^{+}\|_{P}^{2}+2\lambda
L_{\phi}L_{f}\delta\|\varphi_{k}^{+}\|+\lambda(L_{\phi}L_{f}\delta)^{2},$ (19)

where $\|x-x_{k}\|\leq\delta$, $L_{\phi}$ and $L_{f}$ are the Lipschitz
constants of the dictionary functions $\phi$ and dynamics $f$, respectively,
and $\lambda$ is greater than the largest eigenvalue of $P$ supported by
(18c). Similarly, the Lyapunov-like function (9) can be lower-bounded at any
state $x$ with respect to the adjacent sample $x_{k}$ by

$\|\varphi\|_{P}^{2}\geq\|\varphi_{k}\|_{P}^{2}-2\|\varphi-\varphi_{k}\|_{P}\|\varphi_{k}\|_{P}.$

Again, we use Lipschitz continuity to show

$\|\varphi\|_{P}^{2}\geq\|\varphi_{k}\|_{P}^{2}-2\lambda
L_{\phi}\delta\|\varphi_{k}\|.$ (20)

Therefore, the Lyapunov-like function (9) satisfies (8b), where $\beta=0$
i.e.,

$\|\varphi^{+}\|_{P}^{2}-(1-\gamma)\|\varphi\|_{P}^{2}\leq\gamma,$

at any state $x\in\mathcal{X}$ since (19) and (20) hold, and (18b) is enforced
where $\epsilon_{k}$ satisfies (16). Thus, $\|\varphi\|_{P}^{2}\leq 1$ follows
$\|\varphi^{+}\|_{P}^{2}\leq 1$ i.e., trajectories starting at any point
$x\in\mathcal{O}({\bar{r}})$ stay in $\mathcal{O}({\bar{r}})$.

Theorem 1 shows that the computed pi set $\mathcal{O}({\bar{r}})$ is a sub-
maximal pi set, since $\mathcal{O}({\bar{r}})\subseteq\mathcal{X}$ follows
$\mathcal{O}({\bar{r}})\subseteq\mathcal{O}_{\infty}({\bar{r}})$ by
definition, and the cost (18a) rewards enlarging the level-sets of the
Lyapunov-like function $V$ (9) as it minimizes a measure on $V$.

Although Theorem 1 guarantees that the solution to (18) produces a pi set, it
is not guaranteed that (18) will have a solution. The feasibility of sdp (18)
depends on several factors including the existence of the Lyapunov-like
function in the first place, the choice of the lift functions $\phi$, the
data-driven uncertainty compensation $\epsilon_{k}$, and the accuracy of the
equilibrium approximation (15). The following theorem shows that sdp (18) is
feasible if the data (1) is sufficiently dense, the length $T$ of the sampled
state trajectories is sufficiently long, and the lift functions $\phi$ can
describe a Lyapunov-like function of the form (9) that satisfies (8b) with
nonzero $\beta$.

######  Theorem 2.

Let Assumption 1 hold. Let $V(x;\bar{r})$ from (9) satisfy (8b) for some
$0<\beta<1$, where $x_{\infty}(\bar{r})$ is the true equilibrium. There exists
$0<\gamma\leq 1$, $T>0$, and $\delta>0$ such that the Lyapunov-like condition
(18b) holds for the data (1), where $x_{\infty}(\bar{r})$ is approximated by
(15).

Let
$\epsilon_{r}=\frac{1}{N_{T}}\sum_{j=1}^{N_{T}}{\textbf{x}_{j}^{i}}_{T}-x_{\infty}(\bar{r}^{i})$
denote the equilibrium approximation error. We bound the left hand side of the
Lyapunov-like condition (18b) in terms of the exact $x_{\infty}(\bar{r}^{i})$,

$\langle
P^{\prime},\psi_{k}\rangle\leq\|\phi(x_{k}^{+})-\phi(x_{\infty}+\epsilon_{r})\|_{P}^{2}-\|\phi(x_{k})-\phi(x_{\infty}+\epsilon_{r})\|_{P}^{2}.$

By triangle inequality and the Lipschitz continuity of $\phi$ we can write

$\langle P^{\prime},\psi_{k}\rangle\leq
V(x_{k}^{+};\bar{r})-(1-\gamma)V(x_{k};\bar{r})+2\lambda
L_{\phi}\|\epsilon_{r}\|^{2}\leq\gamma-\epsilon_{k}\lambda,$

where, $\lambda$ is the greatest eigenvalue of $P$, and with a slight abuse of
notation, $\psi_{k}$ uses the approximated equilibrium (15) while $V$ uses the
exact equilibrium $x_{\infty}$. The final inequality holds where

$\lambda(2L_{\phi}\|\epsilon_{r}\|^{2}+\epsilon_{k})<\beta.$

This can hold where for instance $\epsilon_{k}<\frac{\beta}{2\lambda}$ by (16)
with some data density $\delta<\\!\\!<\|\varphi_{k}\|$, and the equilibrium
approximation is accurate enough $\|\epsilon_{r}\|^{2}<\frac{\beta}{4\lambda
L_{\phi}}$ by some $T>\\!\\!>0$, since the system is asymptotically iss by
Assumption 1(b).

The exact evaluation of the maximal pi set
$\mathcal{O}({\bar{r}})=\mathcal{O}_{\infty}({\bar{r}})$ depends on the choice
of the dictionary functions $\phi$. Since $\mathcal{O}_{\infty}({\bar{r}})$ is
unknown, an active learning method can be adopted to maximize the volume of
$\mathcal{O}({\bar{r}})$ by changing the dictionary functions $\phi$.

###  3.4 Computational Tractability

Solving sdp (18) is computationally expensive, especially when we use a high-
dimensional lifting $\phi$. In this section, we resolve this issue by relaxing
sdp (18) into a lp by using a series of positive definite matrices
$W_{j}\in\mathbb{R}^{N_{\phi}\times N\phi},j=1,\cdots,N_{W}$, where $0\prec
W_{j}\preceq\frac{1}{M}(\lambda I_{N_{\phi}}-cc^{\top\\!})$, to approximate
the solution to the sdp (18) by

$P=cc^{\top\\!}+\sum_{j=1}^{N_{W}}\alpha_{j}W_{j}.$ (21)

Then, the sdp (18) reduces to the following lp that solves for the vector
$\alpha^{\star}=[\alpha_{1},\cdots,\alpha_{N_{W}}]^{\top\\!}$.

$\displaystyle\alpha^{\star}=\arg\min~{}$ $\displaystyle d^{\top\\!}\alpha$
(22a) s.t. $\displaystyle A\alpha+b\leq 0,$ (22b) $\displaystyle
0\leq\alpha_{j}\leq 1,$ $\displaystyle j=1,\cdots,N_{W},$ (22c)

where the entries of the matrix $A$ and vectors $b$ and $d$ are

$\displaystyle[A]_{k,j}$ $\displaystyle=\langle W_{j},\psi_{k}\rangle$
$\displaystyle[b]_{k}$
$\displaystyle=(c^{\top\\!}\varphi_{k}^{+})^{2}-(c^{\top\\!}\varphi_{k})^{2}-\gamma+\epsilon_{k}\lambda$
$\displaystyle[d]_{j}$ $\displaystyle=\langle W_{j},\Psi\rangle$

respectively. We note that the vector $d$ is defined such that the costs of lp
(22) and sdp (18) be equivalent for the $P$ of the form (21). Moreover, the
following proposition shows that $P$ from (21) is a feasible solution to sdp
(18).

######  Proposition 2.

Let $\alpha$ be the optimal solution of (22). Then, the matrix $P$ (21) is a
feasible solution for the sdp (18),

Condition (22b) enforces $\sum_{j=1}^{N_{W}}[A]_{k,j}\alpha_{j}+[b]_{k}\leq 0$
for $k=1,\cdots,N_{s}$, which is equivalent to (18b) with $P$ from (21).
Condition (22c) enforces $cc^{\top\\!}\preceq P\preceq\lambda I_{N_{\phi}}$
(18c) by triangle inequality where $0\prec W_{j}\preceq\frac{1}{M}(\lambda
I_{N_{\phi}}-cc^{\top\\!})$ in (21). Proposition 2 allows us to avoid the
computational complexity of the sdp (18), which is proportional to the square
of the lift dimension $N_{\phi}^{2}$. Consequently, we can use more nonlinear
functions $\phi$ to capture the nonlinearities of the system in exchange for a
sub-optimal, but feasible, solution to sdp (18). Since Theorem 1 ensures that
feasible solutions of (18) are pi the feasible relaxation (22) also produces
provable pi sets.

##  4 Numerical Examples

In this section, we present numerical results that validate our data-driven
invariant set synthesis procedure. Furthermore, we empirically demonstrate
that the data-driven rg (4) enforces constraints.

###  4.1 Illustrative Example

To illustrate our data-driven rg, we use the following second-order linear
time-invariant (lti) system

$\dot{x}=\begin{bmatrix}0&1\\\
-\omega^{2}&-2\zeta\omega\end{bmatrix}x+\begin{bmatrix}0\\\
\omega^{2}\end{bmatrix}r,~{}~{}x\in\mathcal{X},~{}r\in\mathcal{R}\\\ $ (23)

where $\mathcal{X}=\\{x|~{}y\\!=\\![1,~{}0]x\in[-1,~{}1]\\}$ is the state
constraint set, $\mathcal{R}=[-1.2,~{}1.2]$ is the reference constraint set,
$\omega=5$ is the natural frequency, and $\zeta=0.1$ is the damping ratio of
the system (23). The system is under-damped, thus the overshoot can cause
constraint violations even when the steady-state $x_{\infty}$ lies in the
constraint set $x_{\infty}\in\mathcal{X}$. We purposely chose a linear system
with polytopic constraints (23) since there exist standard model-based methods
[5] to compute the actual maximal admissible set $\Lambda_{\infty}$ for this
problem to compare with our data-driven admissible set $\bar{\Lambda}$.

_Data generation:_ We discretized the system (23) with sampling time $\Delta
t=0.1s$ to generate the data (1) needed for Algorithm 2. For every reference
set-point $\bar{r}_{i}=-1.2+0.02(i-1)\in[-1.2,~{}1.2]$ where $i\in\mathbb{N}$,
we generated $N_{T}=5$ trajectories of length $T=40s$ with random initial
conditions. For the dictionary functions $\phi$, a “complete” set of thin-plat
radial basis functions are considered, centers of which are located on a
$14\times 14$ rectangular grid of states i.e., position and velocity that
encompass the data samples. We used the data to approximate $P$ by solving lp
(22), where $\gamma=0$.

Fig (2) shows the comparison of true maximal and estimated admissible sets,
$\Lambda_{\infty}$ and $\bar{\Lambda}$ from different angles. For instance,
Fig (2).(a) shows the the maximal ci set $\mathcal{C}$ and our data-driven ci
set (5), which are the projection of the admissible sets onto state plane.

Fig. 3 shows the performance of our data-driven rg by Algorithm 1 for the
output position $y=[1,0]x$ of the system (23), where it uses the admissible
set $\bar{\Lambda}$ to adjust the desired reference to enforce the
constraints.

(a) Top view

(b) Front view

(c) Side view

Figure 2: The blue polyhedron is the maximal admissible set $\Lambda_{\infty}$
of the system (23) and the orange balloon like set is our data driven
admissible set $\bar{\Lambda}$, produced by Algorithm (2). Figures are the
projection of the corresponding sets onto (a) $x_{1},x_{2}$ (b) $x_{1},r$, and
(c) $x_{2},r$ planes, where $x_{1}$, $x_{2}$, and $r$ are respectively
position, velocity and reference. Figure 3: Reference tracking and constraint
enforcement for the linear system (23) by our data-driven rg (Algorithm 1).

###  4.2 Autonomous Lane-Keeping

In this section, we showcase the capability of our data-driven rg in handling
nonlinear systems through an autonomous lane-keeping problem. The following
kinematic bicycle model from [24] is borrowed for the vehicle dynamics

$\displaystyle\dot{y}$ $\displaystyle=v\sin(\theta+\beta),\qquad y\in[-2,2]$
(24a) $\displaystyle\dot{\theta}$ $\displaystyle=\frac{v}{l}\sin(\beta),$
(24b)

where the state $x=[y,\theta]^{\top\\!}$ includes the lateral position $y$ of
the center of mass and the inertial heading angle $\theta$, speed $v$ is the
nominal speed of the vehicle, and $l=1.6m$ is the distance from the center of
mass of the vehicle to the rear axle. To facilitate plotting the admissible
set, we will consider the angle of the velocity with respect to the
longitudinal axis of the car $\beta$ as the control input. Thus, we can plot
the $2$ states $y$ and $\theta$ and reference $r$ on a single plot. The state
constraint set $\mathcal{X}=\\{x|~{}[1,0]x\in[-2,2]\\}$ are the allowable
maneuver of the vehicle from the center of the two-lane road, where the widths
of the vehicle and each lane are $1.6m$ and $3.6m$, respectively.

Considering $\beta$ as the control input, we can design a locally stabilizing
controller by linearizing (24) as follows

$\dot{x}=\begin{bmatrix}0&v\\\ 0&0\end{bmatrix}x+\begin{bmatrix}v\\\
v/l\end{bmatrix}\beta.$ (25)

We used linear–quadratic regulator (lqr) with the weights
$Q=\text{diag}([100,10])$ and $R=1$ for state tracking $x\to[r,0]^{\top\\!}$
at the nominal speed $v=20~{}m/s\approx 45mph$. Although the closed-loop
system is almost over-damped under the lqr controller in the nominal speed, it
shows under-damped behavior at higher speed. Our goal is to guarantee
constraint enforcement both at the nominal and higher speed $v=27~{}m/s\approx
60mph$ using the lqr inner-controller designed for the nominal speed.

_Data generation:_ For each cases of nominal and higher speed, we simulated
the kinetic bicycle model (24) using ode45 solver in matlab during the
sampling period $\Delta t=0.1s$ to generate the data needed for Algorithm 2.
For every reference set-point
$\bar{r}_{i}=-1.2+0.02(i-1)\bar{r}\in[-2.4,~{}2.4]$ where $i\in\mathbb{N}$, we
generated $N_{T}=5$ trajectories of length $T=4s$ with random initial
conditions with the nominal lqr inner-controller. For the dictionary functions
$\phi$, similar to the previous example, a complete set of thin-plat radial
basis functions are considered of which centers are located on a $14\times 14$
rectangular grid that encompass the data samples. We used the data to
approximate $P$ by solving lp (22), where $\gamma=0$.

Fig. 4 shows the admissible set $\bar{\Lambda}$ generated by our data-driven
method. Each slice shows the pi set $\mathcal{O}(\bar{r}_{i})$ associated with
a particular reference $\bar{r}_{i}$. Since the system is nonlinear, the true
maximal ci set $\Lambda_{\infty}$ is not known.

Figure 4: Admissible set for the nonlinear kinetic bicycle model (24)
generated by Algorithm 2. Each slice is a pi set $\mathcal{O}(\bar{r})$ for
the constant $\bar{r}$. Figure 5: Constraint enforcement in tracking of
infeasible references for the nonlinear kinetic bicycle model (24) at nominal
speed $v=45mph$. Figure 6: Constraint enforcement for the nonlinear kinetic
bicycle model (24) in an overtaking maneuver. The speed changes from $45mph$
to $60mph$ at $t=4s$

We conducted two experiments for this problem. In the first experiment the
speed is set to the nominal value $v=20m/s$. The desired reference $r_{d}$
commands the position of the center of the car too close to the edge of the
road which will drive the car out of the road. Fig. 5 shows the capability of
our data-driven rg in handling these infeasible reference commands.

In the second experiment, we simulated an overtaking maneuver where the speed
is changed from the nominal $20m/s$ to higher $27m/s$ at $t=4s$. We set the
reference commands inside the constraint set. However, an underdamped response
may still cause constraint violation. Fig. 6 shows our rg can effectively take
the system to the desired position, while enforcing constraints. In the first
lane change, the adjustment on the desired reference $r_{d}$ is minimal, since
the system shows an almost overdamped behavior at nominal speed. Conversely,
rg adjusts the desired reference $r_{d}$ in the second lane change such that
the nonlinear underdamped-like response of the system becomes overdamped,
resulting in constraint satisfaction as well as a smoother ride.

##  5 Conclusion

Our data-driven method creates pi sets directly from data, which is attractive
when the system of interest is not mathematically accessible. We also showed
that with proper setting, these pi sets can be used to create ci and
admissible sets that help solving control problems such as rgs. The advantages
of our pi set computation approach include that it can be offline, and it can
be reduced to a lp. The pi sets have an explicit lifted-quadratic form
suitable for design in other applications, requiring only a mild condition of
Lipschitz continuity on the plant and lifting functions. However, we are
interested to further relax the Lipschitz continuity condition in a
probabilistic paradigm in future works. Although this can only provide
probabilistic guarantee of correctness, it can cover more general classes of
nonlinear systems, make possible the use of broader tools for lifting
functions, and requires fewer data samples.

## References

  * [1] Y. Zheng, B. Shyrokau, and T. Keviczky, “3dop: Comfort-oriented motion planning for automated vehicles with active suspensions,” in _2022 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2022, pp. 390–395.
  * [2] E. Garone, S. Di Cairano, and I. Kolmanovsky, “Reference and command governors for systems with constraints: A survey on theory and applications,” _Automatica_ , vol. 75, pp. 306–328, 2017.
  * [3] D. Masti, V. Breschi, S. Formentin, and A. Bemporad, “Direct data-driven design of neural reference governors,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 4955–4960.
  * [4] S. Nakano, T. W. Nguyen, E. Garone, T. Ibuki, and M. Sampei, “Explicit reference governor on so (3) for torque and pointing constraint management,” _Automatica_ , vol. 155, p. 111103, 2023.
  * [5] D. Bertsekas, “Infinite time reachability of state-space regions by using feedback control,” _IEEE Transactions on Automatic Control_ , vol. 17, no. 5, pp. 604–613, 1972.
  * [6] E. G. Gilbert and K. T. Tan, “Linear systems with state and control constraints: The theory and application of maximal output admissible sets,” _IEEE Transactions on Automatic control_ , vol. 36, no. 9, pp. 1008–1020, 1991.
  * [7] P. Wang, W. Zhang, X. Zhang, and S. S. Ge, “Robust invariance-based explicit reference control for constrained linear systems,” _Automatica_ , vol. 143, p. 110433, 2022.
  * [8] L. Burlion, R. Schieni, and I. V. Kolmanovsky, “A reference governor for linear systems with polynomial constraints,” _Automatica_ , vol. 142, p. 110313, 2022.
  * [9] S. K. Mulagaleti, A. Bemporad, and M. Zanon, “Data-driven synthesis of robust invariant sets and controllers,” _IEEE Control Systems Letters_ , vol. 6, pp. 1676–1681, 2021.
  * [10] E. Gilbert and I. Kolmanovsky, “Nonlinear tracking control in the presence of state and control constraints: a generalized reference governor,” _Automatica_ , vol. 38, no. 12, pp. 2063–2073, 2002.
  * [11] M. Korda, “Computing controlled invariant sets from data using convex optimization,” _SIAM Journal on Control and Optimization_ , vol. 58, no. 5, pp. 2871–2899, 2020.
  * [12] A. Abate, D. Ahmed, M. Giacobbe, and A. Peruffo, “Formal synthesis of lyapunov neural networks,” _IEEE Control Systems Letters_ , vol. 5, no. 3, pp. 773–778, 2020.
  * [13] S. A. Deka and D. V. Dimarogonas, “Supervised learning of lyapunov functions using laplace averages of approximate koopman eigenfunctions,” _IEEE Control Systems Letters_ , 2023.
  * [14] Y.-C. Chang, N. Roohi, and S. Gao, “Neural lyapunov control,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [15] A. Ribeiro, A. Fioravanti, A. Moutinho, and E. de Paiva, “Nonlinear state-feedback design for vehicle lateral control using sum-of-squares programming,” _Vehicle System Dynamics_ , vol. 60, no. 3, pp. 743–769, 2022.
  * [16] S. K. Mulagaleti, A. Bemporad, and M. Zanon, “Data-driven synthesis of robust invariant sets and controllers,” _IEEE Control Systems Letters_ , vol. 6, pp. 1676–1681, 2021.
  * [17] A. Cotorruelo, M. Hosseinzadeh, D. R. Ramirez, D. Limon, and E. Garone, “Reference dependent invariant sets: Sum of squares based computation and applications in constrained control,” _Automatica_ , vol. 129, p. 109614, 2021.
  * [18] S. M. Richards, F. Berkenkamp, and A. Krause, “The lyapunov neural network: Adaptive stability certification for safe learning of dynamical systems,” in _Conference on Robot Learning_. PMLR, 2018, pp. 466–476.
  * [19] J. Osorio, M. Santillo, J. B. Seeds, M. Jankovic, and H. R. Ossareh, “A novel reference governor approach for constraint management of nonlinear systems,” _Automatica_ , vol. 146, p. 110554, 2022.
  * [20] S. Kurt and L. Martin, “Flatness of discrete-time systems, a simple approach,” _arXiv preprint arXiv:2303.05158_ , 2023.
  * [21] F. B. Mathiesen, S. C. Calvert, and L. Laurenti, “Safety certification for stochastic systems via neural barrier functions,” _IEEE Control Systems Letters_ , vol. 7, pp. 973–978, 2022.
  * [22] J. Anderson and A. Papachristodoulou, “Advances in computational lyapunov analysis using sum-of-squares programming.” _Discrete & Continuous Dynamical Systems-Series B_, vol. 20, no. 8, 2015.
  * [23] A. Chakrabarty, D. K. Jha, G. T. Buzzard, Y. Wang, and K. G. Vamvoudakis, “Safe approximate dynamic programming via kernelized lipschitz estimation,” _IEEE transactions on neural networks and learning systems_ , vol. 32, no. 1, pp. 405–419, 2020.
  * [24] R. Rajamani, _Vehicle dynamics and control_. Springer Science & Business Media, 2011.

  *[mpc]: model predictive control
  *[rg]: reference governor
  *[rgs]: reference governor
  *[ci]: controlled invariant
  *[pi]: positively invariant
  *[sdp]: semi-definite program
  *[sos]: sum of squares
  *[iss]: input-to-state stable
  *[ip]: integer program
  *[lp]: linear program
  *[lti]: linear time-invariant
  *[lqr]: linear–quadratic regulator