petitRADTRANS (Sections 6 and 7) and developed a new retrieval approach that
can lead to a robust characterization for cloudy self-luminous atmospheres
(Section 8.3). Specifically, we incorporate the radiative-convective
equilibrium temperature profiles as parameterized priors on the planet’s
thermal structure during the retrievals. This novel approach is enabled by
constraining the temperature-pressure profiles via the temperature gradient
($d\ln{T}/d\ln{P}$), a departure from previous studies that solely modeled the
temperature.

  5. $\bullet$

We have analyzed all published emission spectra and photometry of AF Lep b,
spanning $0.9-4.2$ $\mu$m. Multiple runs are conducted using different
portions of the planet’s spectrophotometry, along with different priors on the
planet’s mass and radius. These retrievals consistently suggest that AF Lep b
likely has a metal-enriched atmosphere ([Fe/H]$>1.0$ dex or $Z_{\rm
planet}/Z_{\rm star}=75^{+94}_{-42}$) compared to the metallicity of its host
star. Our analysis also determines $T_{\rm eff}\approx 800$ K,
$\log{(g)}\approx 3.7$ dex, and the presence of silicate clouds and the dis-
equilibrium chemistry in the atmosphere of AF Lep b.

  6. $\bullet$

The potential metal enrichment of AF Lep b might be explained by planetesimal
accretion, giant impacts, and/or core erosion. The first process also
coincides with the presence of a debris disk in the system located at $40-60$
au. At its observed orbit, AF Lep b has a sufficient mass to dynamically
excite the debris disk, causing planetesimal scattering and bombardment onto
the planet’s atmosphere (Section 8.1).

  7. $\bullet$

If AF Lep b formed a few Myr later than its host star, then the planet’s
retrieved effective temperature and bolometric luminosity would be consistent
with predictions of hot-start evolution models. When compared to predictions
of cold-start evolution models, our retrieved $T_{\rm eff}$ and $L_{\rm bol}$
appear to be much higher regardless of whether the planet is coeval with the
host star or formed 10 Myr later (Section 8.1).

  8. $\bullet$

If giant impacts and planetary mergers occurred during the evolution history
of AF Lep b, these processes might act as a mechanism of rejuvenation, by
altering the planet’s entropy state. In this scenario, the noted discrepancies
bewteen retrieved and evolution-based $T_{\rm eff}$ and $L_{\rm bol}$ of AF
Lep b can be potentially resolved.

  9. $\bullet$

We have refined the spectral type of AF Lep b to be L8–T3. Compared to higher-
mass brown dwarfs with similar spectral types, AF Lep b has $450-600$ K colder
$T_{\rm eff}$ and $0.6-0.9$ dex fainter $L_{\rm bol}$. These peculiar
properties of AF Lep b are likely linked to its very low surface gravity and
high atmospheric metallicity. Notably, AF Lep b is the coldest object with
suggested evidence of silicate clouds to date (Section 8.2).

  10. $\bullet$

Our analysis of AF Lep b also leads to several valuable lessons that may have
broader implications for the atmospheric study of directly imaged exoplanets.
Specifically, we highlight the values of (1) performing the end-to-end
reduction for the same set of direct imaging data by using different
pipelines, (2) measuring the planet’s dynamical mass based on orbital
monitoring, (3) contextualizing the planets’ atmospheric properties by using
evolution models, (4) acquiring the planet’s spectrophotometry with a broad
wavelength coverage, and (5) studying the planet’s top-of-atmosphere
inhomogeneities based on the variability monitoring.

For future work, it is essential to acquire spectrophotometry of AF Lep b with
a higher S/N and/or spectral resolution over a broader wavelength range (e.g.,
VLTI/GRAVITY). Also, consistent reduction procedures should be applied to the
existing SPHERE data (from Mesa et al. 2023 and De Rosa et al. 2023) to
improve the quality of the emission spectra. These efforts will help improve
the accuracy and precision of the inferred atmospheric properties and of AF
Lep b, allowing for a reassessment of the formation pathway of this remarkable
planetary system.

Z. Z. and P. M. thank Michiel Min for the idea of parameterizing the thermal
structure using the temperature gradient. Z. Z. thanks Artem Aguichine,
William Balmer, Bertram Bitsch, Douglas Lin, Evert Nasedkin, Kazumasa Ohno,
and Daniel Thorngren for very helpful discussions. We thank the referee for
suggestions that improved the manuscript. Support for this work was provided
by NASA through the NASA Hubble Fellowship grant HST-HF2-51522.001-A awarded
by the Space Telescope Science Institute, which is operated by the Association
of Universities for Research in Astronomy, Inc., for NASA, under contract
NAS5-26555. B.P.B. acknowledges support from the National Science Foundation
grant AST-1909209, NASA Exoplanet Research Program grant 20-XRP20$\\_$2-0119,
and the Alfred P. Sloan Foundation. This work has benefited from The
UltracoolSheet at http://bit.ly/UltracoolSheet, maintained by Will Best, Trent
Dupuy, Michael Liu, Rob Siverd, and Zhoujian Zhang, and developed from
compilations by Dupuy & Liu (2012), Dupuy & Kraus (2013), Liu et al. (2016),
Best et al. (2018), and Best et al. (2021). This work has made use of data
from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular, the institutions
participating in the Gaia Multilateral Agreement.

##  Appendix A Literature Comparison for our Stellar parameters of AF Lep A

As shown in Figure 16, our inferred mass, iron abundance, effective
temperature, surface gravity, and radius of AF Lep A (Table 2.2) are generally
consistent with those determined in the literature. When collecting the
literature values, we exclude the studies that cataloged or averaged the
stellar parameters measured by other work.

In this work, $M_{A}$ is a key parameter for constraining the orbit and
dynamical mass of AF Lep b. Several studies measured slightly higher $M_{A}$
than our adopted value, including Valenti & Fischer (2005), Takeda et al.
(2007), Wright et al. (2011), Kervella et al. (2019), and Pearce et al.
(2022); the latter four studies did not measure the stellar [Fe/H]. As
detailed below, the different $M_{A}$ estimates between their and our works
are likely due to the systematic differences of the adopted stellar model
spectra and isochrones (with different atomic/molecular line lists), input
photometry (single band versus multiple bands), the assumption about the
stellar age, and whether all input parameters are simultaneously or separately
constrained by the stellar models.

$\bullet$ Valenti & Fischer (2005) estimated the bolometric luminosity of AF
Lep A via $V$-band magnitude with a bolometric correction and derived other
stellar parameters by fitting the observed spectrum using the Kurucz (1992)
model atmospheres; they then updated the stellar mass by modeling the $L_{\rm
bol}$ and $T_{\rm eff}$ via the Yonsei-Yale isochrones (Demarque et al.,
2004). Their resulting $M_{A}$, [Fe/H]A, $T_{\rm eff,A}$, and $\log{(g_{A})}$
are all among the highest (and their derived $R_{A}$ is among the lowest) in
the literature.

$\bullet$ The analysis of Takeda et al. (2007) is tied to that of Valenti &
Fischer (2005). Briefly, Takeda et al. (2007) derived the mass and radius of
AF Lep A (along with about $1000$ cool stars) by modeling this object’s
parallax, $V$-band magnitude, and the Valenti & Fischer (2005) $T_{\rm eff}$,
$\log{(g)}$, [Fe/H] measurements using the Yale Rotational Evolution Code.

$\bullet$ Wright et al. (2011) first estimated a bolometric luminosity based
on the $V$-band magnitude and a bolometric correction guided by the object’s
$V-K_{s}$ color. Assuming AF Lep A has a much older age of $1$ Gyr (compared
to its age of $24\pm 3$ Myr), they then used the bolometric luminosity and the
Siess et al. (2000) isochrones to derive stellar $M_{A}$, $R_{A}$, and $T_{\rm
eff,A}$.

$\bullet$ Kervella et al. (2019) estimated the radius of AF Lep A using the
$V$-band magnitude, $V-K_{s}$ color, and the parallax. Then they derived
$M_{A}$ by using the Girardi et al. (2000) isochrones.

$\bullet$ Pearce et al. (2022) determined the effective temperature and
bolometric luminosity of AF Lep A by fitting the PHOENIX (Husser et al., 2013)
model atmospheres (with a blackbody component to account for the debris disk)
to the observed multi-band photometry. They then modeled these inferred
$T_{\rm eff,A}$ and $L_{\rm bol,A}$ using the MIST evolution models to derive
$M_{A}$.

##  Appendix B Parameter Posteriors and Data–Model Comparison of Our Orbit
Analysis

In Figure 17, we present the posteriors for masses of AF Lep A and b, as well
as the planet’s orbital parameters. Figure 18 compares our input observational
data with the fitted orbits.

##  Appendix C Evolution-based Properties of AF Lep b Assuming a Younger Age
than the Host Star

Table 7 presents the physical properties of AF Lep b inferred by various
evolution models. Input parameters of this analysis include the directly
measured mass ($M=2.8^{+0.6}_{-0.5}$ MJup; Section 4) and an age of $14\pm 3$
My that is 10 Myr younger than that of its host star, assuming a slightly
later epoch of planet formation.

##  Appendix D Retrievals on All $K1/K2/L^{\prime}$ Photometry with Different
Spectroscopic Datasets and Parameter Priors

Figure 19 presents the results of all retrieval runs that incorporate all the
$K1/K2/L^{\prime}$ photometry, with different sets of spectra and parameter
priors. As described in Section 7, the three spectral sets include: (1) the
Mesa et al. (2023) spectrum, (2) the De Rosa et al. (2023) spectrum, and (3)
both spectra. The three sets of parameter priors include: (1) the default
priors summarized in Table 5, (2) a constrained prior on the planet’s mass and
default priors for the remaining parameters, and (3) constrained priors on
both the planet’s mass and radius with default priors for the remaining
parameters. Retrievals with constrained priors on both $M$ and $R$ are already
shown in Figures 7–9, so Figure 19 presents the results for the remaining six
retrieval runs. Median values and confidence intervals of all parameters
inferred by these retrievals are summarized in Table 6.

##  Appendix E Retrievals on Different Spectroscopic Datasets and Parameter
Priors of AF Lep b, with Photometry Excluded

Figure 20 presents the results of the additional nine retrieval runs that
incorporate different sets of spectra and parameter priors (same as described
in Appendix D), but with all the $K1/K2/L^{\prime}$ photometry of AF Lep b
excluded.

## References

  * Ackerman & Marley (2001) Ackerman, A. S., & Marley, M. S. 2001, ApJ, 556, 872, doi: 10.1086/321540
  * Ahrer et al. (2023) Ahrer, E.-M., Stevenson, K. B., Mansfield, M., et al. 2023, Nature, 614, 653, doi: 10.1038/s41586-022-05590-4
  * Alexander et al. (2014) Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 475–496, doi: 10.2458/azu_uapress_9780816531240-ch021
  * Ali-Dib et al. (2022) Ali-Dib, M., Cumming, A., & Lin, D. N. C. 2022, MNRAS, 509, 1413, doi: 10.1093/mnras/stab3008
  * Alibert et al. (2005) Alibert, Y., Mousis, O., & Benz, W. 2005, ApJ, 622, L145, doi: 10.1086/429662
  * Allard et al. (2001) Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357, doi: 10.1086/321547
  * Allard et al. (2012) Allard, F., Homeier, D., & Freytag, B. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2765, doi: 10.1098/rsta.2011.0269
  * Allard et al. (2016) Allard, N. F., Spiegelman, F., & Kielkopf, J. F. 2016, A&A, 589, A21, doi: 10.1051/0004-6361/201628270
  * Allard et al. (2019) Allard, N. F., Spiegelman, F., Leininger, T., & Molliere, P. 2019, A&A, 628, A120, doi: 10.1051/0004-6361/201935593
  * Allende Prieto & Lambert (1999) Allende Prieto, C., & Lambert, D. L. 1999, A&A, 352, 555, doi: 10.48550/arXiv.astro-ph/9911002
  * Ammons et al. (2006) Ammons, S. M., Robinson, S. E., Strader, J., et al. 2006, ApJ, 638, 1004, doi: 10.1086/498490
  * Anders et al. (2019) Anders, F., Khalatyan, A., Chiappini, C., et al. 2019, A&A, 628, A94, doi: 10.1051/0004-6361/201935765
  * Anders et al. (2022) Anders, F., Khalatyan, A., Queiroz, A. B. A., et al. 2022, A&A, 658, A91, doi: 10.1051/0004-6361/202142369
  * Anderson & Francis (2012) Anderson, E., & Francis, C. 2012, Astronomy Letters, 38, 331, doi: 10.1134/S1063773712050015
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Atreya et al. (2016) Atreya, S. K., Crida, A., Guillot, T., et al. 2016, arXiv e-prints, arXiv:1606.04510, doi: 10.48550/arXiv.1606.04510
  * August et al. (2023) August, P. C., Bean, J. L., Zhang, M., et al. 2023, arXiv e-prints, arXiv:2305.07753, doi: 10.48550/arXiv.2305.07753
  * Azzam et al. (2016) Azzam, A. A. A., Tennyson, J., Yurchenko, S. N., & Naumenko, O. V. 2016, MNRAS, 460, 4063, doi: 10.1093/mnras/stw1133
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147, doi: 10.3847/1538-3881/abd806
  * Baraffe et al. (2002) Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 2002, A&A, 382, 563, doi: 10.1051/0004-6361:20011638
  * Baraffe et al. (2003) Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701, doi: 10.1051/0004-6361:20030252
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42, doi: 10.1051/0004-6361/201425481
  * Beatty et al. (2018) Beatty, T. G., Morley, C. V., Curtis, J. L., et al. 2018, AJ, 156, 168, doi: 10.3847/1538-3881/aad697
  * Bell et al. (2015) Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593, doi: 10.1093/mnras/stv1981
  * Best et al. (2020b) Best, W. M. J., Dupuy, T. J., Liu, M. C., Siverd, R. J., & Zhang, Z. 2020b, The UltracoolSheet: Photometry, Astrometry, Spectroscopy, and Multiplicity for 3000+ Ultracool Dwarfs and Imaged Exoplanets, Zenodo, doi: 10.5281/zenodo.4169085
  * Best et al. (2021) Best, W. M. J., Liu, M. C., Magnier, E. A., & Dupuy, T. J. 2021, AJ, 161, 42, doi: 10.3847/1538-3881/abc893
  * Best et al. (2018) Best, W. M. J., Magnier, E. A., Liu, M. C., et al. 2018, ApJS, 234, 1, doi: 10.3847/1538-4365/aa9982
  * Biersteker & Schlichting (2019) Biersteker, J. B., & Schlichting, H. E. 2019, MNRAS, 485, 4454, doi: 10.1093/mnras/stz738
  * Blanco-Cuaresma (2019) Blanco-Cuaresma, S. 2019, MNRAS, 486, 2075, doi: 10.1093/mnras/stz549
  * Blanco-Cuaresma et al. (2014) Blanco-Cuaresma, S., Soubiran, C., Heiter, U., & Jofré, P. 2014, A&A, 569, A111, doi: 10.1051/0004-6361/201423945
  * Boley & Durisen (2010) Boley, A. C., & Durisen, R. H. 2010, ApJ, 724, 618, doi: 10.1088/0004-637X/724/1/618
  * Bonavita et al. (2022) Bonavita, M., Fontanive, C., Gratton, R., et al. 2022, MNRAS, 513, 5588, doi: 10.1093/mnras/stac1250
  * Borysow (2002) Borysow, A. 2002, A&A, 390, 779, doi: 10.1051/0004-6361:20020555
  * Borysow & Frommhold (1989) Borysow, A., & Frommhold, L. 1989, ApJ, 341, 549, doi: 10.1086/167515
  * Borysow et al. (1989) Borysow, A., Frommhold, L., & Moraldi, M. 1989, ApJ, 336, 495, doi: 10.1086/167027
  * Borysow et al. (2001) Borysow, A., Jorgensen, U. G., & Fu, Y. 2001, J. Quant. Spec. Radiat. Transf., 68, 235, doi: 10.1016/S0022-4073(00)00023-6
  * Borysow et al. (1988) Borysow, J., Frommhold, L., & Birnbaum, G. 1988, ApJ, 326, 509, doi: 10.1086/166112
  * Boucher et al. (2023) Boucher, A., Lafreniére, D., Pelletier, S., et al. 2023, MNRAS, 522, 5062, doi: 10.1093/mnras/stad1247
  * Bowler et al. (2021) Bowler, B. P., Endl, M., Cochran, W. D., et al. 2021, ApJ, 913, L26, doi: 10.3847/2041-8213/abfec8
  * Bowler et al. (2023) Bowler, B. P., Tran, Q. H., Zhang, Z., et al. 2023, AJ, 165, 164, doi: 10.3847/1538-3881/acbd34
  * Brandt (2018) Brandt, T. D. 2018, ApJS, 239, 31, doi: 10.3847/1538-4365/aaec06
  * Brandt (2021) —. 2021, ApJS, 254, 42, doi: 10.3847/1538-4365/abf93c
  * Brandt et al. (2020) Brandt, T. D., Dupuy, T. J., Bowler, B. P., et al. 2020, AJ, 160, 196, doi: 10.3847/1538-3881/abb45e
  * Brandt et al. (2021) Brandt, T. D., Dupuy, T. J., Li, Y., et al. 2021, AJ, 162, 186, doi: 10.3847/1538-3881/ac042e
  * Brogi et al. (2023) Brogi, M., Emeka-Okafor, V., Line, M. R., et al. 2023, AJ, 165, 91, doi: 10.3847/1538-3881/acaf5c
  * Brown-Sevilla et al. (2022) Brown-Sevilla, S. B., Maire, A. L., Mollière, P., et al. 2022, arXiv e-prints, arXiv:2211.14330, doi: 10.48550/arXiv.2211.14330
  * Brown-Sevilla et al. (2023) —. 2023, A&A, 673, A98, doi: 10.1051/0004-6361/202244826
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125, doi: 10.1051/0004-6361/201322971
  * Burgasser et al. (2006) Burgasser, A. J., Burrows, A., & Kirkpatrick, J. D. 2006, ApJ, 639, 1095, doi: 10.1086/499344
  * Burningham et al. (2017) Burningham, B., Marley, M. S., Line, M. R., et al. 2017, MNRAS, 470, 1177, doi: 10.1093/mnras/stx1246
  * Burningham et al. (2009) Burningham, B., Pinfield, D. J., Leggett, S. K., et al. 2009, MNRAS, 395, 1237, doi: 10.1111/j.1365-2966.2009.14620.x
  * Burningham et al. (2021) Burningham, B., Faherty, J. K., Gonzales, E. C., et al. 2021, MNRAS, 506, 1944, doi: 10.1093/mnras/stab1361
  * Burrows et al. (1997) Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856, doi: 10.1086/305002
  * Butler et al. (2017) Butler, R. P., Vogt, S. S., Laughlin, G., et al. 2017, AJ, 153, 208, doi: 10.3847/1538-3881/aa66ca
  * Casagrande et al. (2011) Casagrande, L., Schönrich, R., Asplund, M., et al. 2011, A&A, 530, A138, doi: 10.1051/0004-6361/201016276
  * Chabrier et al. (2000) Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ApJ, 542, 464, doi: 10.1086/309513
  * Chan & Dalgarno (1965) Chan, Y. M., & Dalgarno, A. 1965, Proceedings of the Physical Society, 85, 227, doi: 10.1088/0370-1328/85/2/304
  * Chandler et al. (2016) Chandler, C. O., McDonald, I., & Kane, S. R. 2016, AJ, 151, 59, doi: 10.3847/0004-6256/151/3/59
  * Changeat et al. (2022) Changeat, Q., Edwards, B., Al-Refaie, A. F., et al. 2022, ApJS, 260, 3, doi: 10.3847/1538-4365/ac5cc2
  * Charnay et al. (2018) Charnay, B., Bézard, B., Baudino, J. L., et al. 2018, ApJ, 854, 172, doi: 10.3847/1538-4357/aaac7d
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Claudi et al. (2008) Claudi, R. U., Turatto, M., Gratton, R. G., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II, ed. I. S. McLean & M. M. Casali, 70143E, doi: 10.1117/12.788366
  * Coles et al. (2019) Coles, P. A., Yurchenko, S. N., & Tennyson, J. 2019, MNRAS, 490, 4638, doi: 10.1093/mnras/stz2778
  * Cridland et al. (2020) Cridland, A. J., van Dishoeck, E. F., Alessi, M., & Pudritz, R. E. 2020, A&A, 642, A229, doi: 10.1051/0004-6361/202038767
  * Currie et al. (2023) Currie, T., Brandt, G. M., Brandt, T. D., et al. 2023, Science, 380, 198, doi: 10.1126/science.abo6192
  * Cushing et al. (2006) Cushing, M. C., Roellig, T. L., Marley, M. S., et al. 2006, ApJ, 648, 614, doi: 10.1086/505637
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, II/246
  * Cutri et al. (2021) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR Online Data Catalog, II/328
  * Dahn et al. (2002) Dahn, C. C., Harris, H. C., Vrba, F. J., et al. 2002, AJ, 124, 1170, doi: 10.1086/341646
  * Dalgarno & Williams (1962) Dalgarno, A., & Williams, D. A. 1962, ApJ, 136, 690, doi: 10.1086/147428
  * De Rosa et al. (2023) De Rosa, R. J., Nielsen, E. L., Wahhaj, Z., et al. 2023, A&A, 672, A94, doi: 10.1051/0004-6361/202345877
  * Debras & Chabrier (2019) Debras, F., & Chabrier, G. 2019, ApJ, 872, 100, doi: 10.3847/1538-4357/aaff65
  * Delorme et al. (2017) Delorme, P., Meunier, N., Albert, D., et al. 2017, in SF2A-2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. C. Reylé, P. Di Matteo, F. Herpin, E. Lagadec, A. Lançon, Z. Meliani, & F. Royer, Di, doi: 10.48550/arXiv.1712.06948
  * Demarque et al. (2004) Demarque, P., Woo, J.-H., Kim, Y.-C., & Yi, S. K. 2004, ApJS, 155, 667, doi: 10.1086/424966
  * Do Ó et al. (2023) Do Ó, C. R., O’Neil, K. K., Konopacky, Q. M., et al. 2023, AJ, 166, 48, doi: 10.3847/1538-3881/acdc9a
  * Dohlen et al. (2008) Dohlen, K., Langlois, M., Saisse, M., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II, ed. I. S. McLean & M. M. Casali, 70143L, doi: 10.1117/12.789786
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8, doi: 10.3847/0067-0049/222/1/8
  * Drazkowska et al. (2022) Drazkowska, J., Bitsch, B., Lambrechts, M., et al. 2022, arXiv e-prints, arXiv:2203.09759, doi: 10.48550/arXiv.2203.09759
  * Dupuy & Kraus (2013) Dupuy, T. J., & Kraus, A. L. 2013, Science, 341, 1492, doi: 10.1126/science.1241917
  * Dupuy & Liu (2012) Dupuy, T. J., & Liu, M. C. 2012, ApJS, 201, 19, doi: 10.1088/0067-0049/201/2/19
  * Dupuy et al. (2009) Dupuy, T. J., Liu, M. C., & Ireland, M. J. 2009, ApJ, 692, 729, doi: 10.1088/0004-637X/692/1/729
  * El-Badry et al. (2021) El-Badry, K., Rix, H.-W., & Heintz, T. M. 2021, MNRAS, 506, 2269, doi: 10.1093/mnras/stab323
  * Emsenhuber et al. (2021) Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021, A&A, 656, A69, doi: 10.1051/0004-6361/202038553
  * Fabricius et al. (2021) Fabricius, C., Luri, X., Arenou, F., et al. 2021, A&A, 649, A5, doi: 10.1051/0004-6361/202039834
  * Feltzing et al. (2001) Feltzing, S., Holmberg, J., & Hurley, J. R. 2001, A&A, 377, 911, doi: 10.1051/0004-6361:20011119
  * Feroz & Hobson (2008) Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449, doi: 10.1111/j.1365-2966.2007.12353.x
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601, doi: 10.1111/j.1365-2966.2009.14548.x
  * Feroz et al. (2019) Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, The Open Journal of Astrophysics, 2, 10, doi: 10.21105/astro.1306.2144
  * Filippazzo et al. (2015) Filippazzo, J. C., Rice, E. L., Faherty, J., et al. 2015, ApJ, 810, 158, doi: 10.1088/0004-637X/810/2/158
  * Finnerty et al. (2023) Finnerty, L., Schofield, T., Sappey, B., et al. 2023, AJ, 166, 31, doi: 10.3847/1538-3881/acda91
  * Fletcher et al. (2009) Fletcher, L. N., Orton, G. S., Teanby, N. A., Irwin, P. G. J., & Bjoraker, G. L. 2009, Icarus, 199, 351, doi: 10.1016/j.icarus.2008.09.019
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fortney et al. (2008) Fortney, J. J., Marley, M. S., Saumon, D., & Lodders, K. 2008, ApJ, 683, 1104, doi: 10.1086/589942
  * Fortney & Nettelmann (2010) Fortney, J. J., & Nettelmann, N. 2010, Space Sci. Rev., 152, 423, doi: 10.1007/s11214-009-9582-x
  * Fouesneau et al. (2022) Fouesneau, M., Andrae, R., Dharmawardena, T., et al. 2022, A&A, 662, A125, doi: 10.1051/0004-6361/202141828
  * Franson et al. (2023a) Franson, K., Bowler, B. P., Bonavita, M., et al. 2023a, AJ, 165, 39, doi: 10.3847/1538-3881/aca408
  * Franson et al. (2023b) Franson, K., Bowler, B. P., Zhou, Y., et al. 2023b, ApJ, 950, L19, doi: 10.3847/2041-8213/acd6f6
  * Frelikh et al. (2019) Frelikh, R., Jang, H., Murray-Clay, R. A., & Petrovich, C. 2019, ApJ, 884, L47, doi: 10.3847/2041-8213/ab4a7b
  * Fu et al. (2022) Fu, G., Sing, D. K., Lothringer, J. D., et al. 2022, ApJ, 925, L3, doi: 10.3847/2041-8213/ac4968
  * Gaarn et al. (2023) Gaarn, J., Burningham, B., Faherty, J. K., et al. 2023, MNRAS, 521, 5761, doi: 10.1093/mnras/stad753
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211, doi: 10.48550/arXiv.2208.00211
  * Gao et al. (2020) Gao, P., Thorngren, D. P., Lee, E. K. H., et al. 2020, Nature Astronomy, 4, 951, doi: 10.1038/s41550-020-1114-3
  * Ginzburg & Chiang (2020) Ginzburg, S., & Chiang, E. 2020, MNRAS, 498, 680, doi: 10.1093/mnras/staa2500
  * Girardi et al. (2000) Girardi, L., Bressan, A., Bertelli, G., & Chiosi, C. 2000, A&AS, 141, 371, doi: 10.1051/aas:2000126
  * Gizis et al. (2012) Gizis, J. E., Faherty, J. K., Liu, M. C., et al. 2012, AJ, 144, 94, doi: 10.1088/0004-6256/144/4/94
  * Glebocki & Gnacinski (2005) Glebocki, R., & Gnacinski, P. 2005, VizieR Online Data Catalog, III/244
  * Golimowski et al. (2004) Golimowski, D. A., Leggett, S. K., Marley, M. S., et al. 2004, AJ, 127, 3516, doi: 10.1086/420709
  * Gonzales et al. (2020) Gonzales, E. C., Burningham, B., Faherty, J. K., et al. 2020, ApJ, 905, 46, doi: 10.3847/1538-4357/abbee2
  * Gonzales et al. (2022) —. 2022, ApJ, 938, 56, doi: 10.3847/1538-4357/ac8f2a
  * Gonzales et al. (2021) —. 2021, ApJ, 923, 19, doi: 10.3847/1538-4357/ac294e
  * GRAVITY Collaboration et al. (2020) GRAVITY Collaboration, Nowak, M., Lacour, S., et al. 2020, A&A, 633, A110, doi: 10.1051/0004-6361/201936898
  * Gray et al. (2006) Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006, AJ, 132, 161, doi: 10.1086/504637
  * Greenbaum et al. (2018) Greenbaum, A. Z., Pueyo, L., Ruffio, J.-B., et al. 2018, AJ, 155, 226, doi: 10.3847/1538-3881/aabcb8
  * Guillot & Gladman (2000) Guillot, T., & Gladman, B. 2000, in Astronomical Society of the Pacific Conference Series, Vol. 219, Disks, Planetesimals, and Planets, ed. G. Garzón, C. Eiroa, D. de Winter, & T. J. Mahoney, 475
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951, doi: 10.1051/0004-6361:200809724
  * Hargreaves et al. (2020) Hargreaves, R. J., Gordon, I. E., Rey, M., et al. 2020, ApJS, 247, 55, doi: 10.3847/1538-4365/ab7a1a
  * Hawkins et al. (2020) Hawkins, K., Lucey, M., & Curtis, J. 2020, MNRAS, 496, 2422, doi: 10.1093/mnras/staa1673
  * Heiter et al. (2021) Heiter, U., Lind, K., Bergemann, M., et al. 2021, A&A, 645, A106, doi: 10.1051/0004-6361/201936291
  * Helled et al. (2011) Helled, R., Anderson, J. D., Podolak, M., & Schubert, G. 2011, ApJ, 726, 15, doi: 10.1088/0004-637X/726/1/15
  * Helled & Bodenheimer (2010) Helled, R., & Bodenheimer, P. 2010, Icarus, 207, 503, doi: 10.1016/j.icarus.2009.11.023
  * Helled & Bodenheimer (2014) —. 2014, ApJ, 789, 69, doi: 10.1088/0004-637X/789/1/69
  * Helled et al. (2006) Helled, R., Podolak, M., & Kovetz, A. 2006, Icarus, 185, 64, doi: 10.1016/j.icarus.2006.06.011
  * Helled & Schubert (2009) Helled, R., & Schubert, G. 2009, ApJ, 697, 1256, doi: 10.1088/0004-637X/697/2/1256
  * Helled et al. (2014) Helled, R., Bodenheimer, P., Podolak, M., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 643–665, doi: 10.2458/azu_uapress_9780816531240-ch028
  * Hoch et al. (2023) Hoch, K. K. W., Konopacky, Q. M., Theissen, C. A., et al. 2023, AJ, 166, 85, doi: 10.3847/1538-3881/ace442
  * Høg et al. (2000) Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27
  * Hood et al. (2023) Hood, C. E., Fortney, J. J., Line, M. R., & Faherty, J. K. 2023, arXiv e-prints, arXiv:2303.04885, doi: 10.48550/arXiv.2303.04885
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6, doi: 10.1051/0004-6361/201219058
  * Järvinen et al. (2015) Järvinen, S. P., Arlt, R., Hackman, T., et al. 2015, A&A, 574, A25, doi: 10.1051/0004-6361/201424229
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python.  http://www.scipy.org/
  * Karalidi et al. (2021) Karalidi, T., Marley, M., Fortney, J. J., et al. 2021, ApJ, 923, 269, doi: 10.3847/1538-4357/ac3140
  * Kervella et al. (2019) Kervella, P., Arenou, F., Mignard, F., & Thévenin, F. 2019, A&A, 623, A72, doi: 10.1051/0004-6361/201834371
  * Kervella et al. (2022) Kervella, P., Arenou, F., & Thévenin, F. 2022, A&A, 657, A7, doi: 10.1051/0004-6361/202142146
  * Kirkpatrick et al. (2021) Kirkpatrick, J. D., Gelino, C. R., Faherty, J. K., et al. 2021, ApJS, 253, 7, doi: 10.3847/1538-4365/abd107
  * Konopacky et al. (2013) Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Science, 339, 1398, doi: 10.1126/science.1232003
  * Kurucz (1992) Kurucz, R. L. 1992, in The Stellar Populations of Galaxies, ed. B. Barbuy & A. Renzini, Vol. 149, 225
  * Kurucz (1993) Kurucz, R. L. 1993, SYNTHE spectrum synthesis programs and line data
  * Kuzuhara et al. (2022) Kuzuhara, M., Currie, T., Takarada, T., et al. 2022, ApJ, 934, L18, doi: 10.3847/2041-8213/ac772f
  * Lacy & Burrows (2023) Lacy, B., & Burrows, A. 2023, arXiv e-prints, arXiv:2303.16295, doi: 10.48550/arXiv.2303.16295
  * Lavie et al. (2017) Lavie, B., Mendonça, J. M., Mordasini, C., et al. 2017, AJ, 154, 91, doi: 10.3847/1538-3881/aa7ed8
  * Leconte (2018) Leconte, J. 2018, ApJ, 853, L30, doi: 10.3847/2041-8213/aaaa61
  * Leggett et al. (2007) Leggett, S. K., Marley, M. S., Freedman, R., et al. 2007, ApJ, 667, 537, doi: 10.1086/519948
  * Li et al. (2010) Li, S. L., Agnor, C. B., & Lin, D. N. C. 2010, ApJ, 720, 1161, doi: 10.1088/0004-637X/720/2/1161
  * Lin & Ida (1997) Lin, D. N. C., & Ida, S. 1997, ApJ, 477, 781, doi: 10.1086/303738
  * Lindegren et al. (2021) Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4, doi: 10.1051/0004-6361/202039653
  * Line et al. (2015) Line, M. R., Teske, J., Burningham, B., Fortney, J. J., & Marley, M. S. 2015, ApJ, 807, 183, doi: 10.1088/0004-637X/807/2/183
  * Line et al. (2017) Line, M. R., Marley, M. S., Liu, M. C., et al. 2017, ApJ, 848, 83, doi: 10.3847/1538-4357/aa7ff0
  * Line et al. (2021) Line, M. R., Brogi, M., Bean, J. L., et al. 2021, Nature, 598, 580, doi: 10.1038/s41586-021-03912-6
  * Liu et al. (2016) Liu, M. C., Dupuy, T. J., & Allers, K. N. 2016, ApJ, 833, 96, doi: 10.3847/1538-4357/833/1/96
  * Liu et al. (2007) Liu, M. C., Leggett, S. K., & Chiu, K. 2007, ApJ, 660, 1507, doi: 10.1086/512662
  * Liu et al. (2013) Liu, M. C., Magnier, E. A., Deacon, N. R., et al. 2013, ApJ, 777, L20, doi: 10.1088/2041-8205/777/2/L20
  * Liu et al. (2015) Liu, S.-F., Agnor, C. B., Lin, D. N. C., & Li, S.-L. 2015, MNRAS, 446, 1685, doi: 10.1093/mnras/stu2205
  * Liu et al. (2019) Liu, S.-F., Hori, Y., Müller, S., et al. 2019, Nature, 572, 355, doi: 10.1038/s41586-019-1470-2
  * Lodders & Fegley (2006) Lodders, K., & Fegley, B., J. 2006, in Astrophysics Update 2, ed. J. W. Mason, 1, doi: 10.1007/3-540-30313-8_1
  * Lothringer et al. (2021) Lothringer, J. D., Rustamkulov, Z., Sing, D. K., et al. 2021, ApJ, 914, 12, doi: 10.3847/1538-4357/abf8a9
  * Luck (2018) Luck, R. E. 2018, AJ, 155, 111, doi: 10.3847/1538-3881/aaa9b5
  * Lueber et al. (2022) Lueber, A., Kitzmann, D., Bowler, B. P., Burgasser, A. J., & Heng, K. 2022, ApJ, 930, 136, doi: 10.3847/1538-4357/ac63b9
  * Lunine et al. (1986) Lunine, J. I., Hubbard, W. B., & Marley, M. S. 1986, ApJ, 310, 238, doi: 10.1086/164678
  * MacDonald & Batalha (2023) MacDonald, R. J., & Batalha, N. E. 2023, Research Notes of the American Astronomical Society, 7, 54, doi: 10.3847/2515-5172/acc46a
  * Madhusudhan et al. (2014) Madhusudhan, N., Amin, M. A., & Kennedy, G. M. 2014, ApJ, 794, L12, doi: 10.1088/2041-8205/794/1/L12
  * Madhusudhan et al. (2017) Madhusudhan, N., Bitsch, B., Johansen, A., & Eriksson, L. 2017, MNRAS, 469, 4102, doi: 10.1093/mnras/stx1139
  * Maíz Apellániz & Weiler (2018) Maíz Apellániz, J., & Weiler, M. 2018, A&A, 619, A180, doi: 10.1051/0004-6361/201834051
  * Mamajek (2009) Mamajek, E. E. 2009, in American Institute of Physics Conference Series, Vol. 1158, Exoplanets and Disks: Their Formation and Diversity, ed. T. Usuda, M. Tamura, & M. Ishii, 3–10, doi: 10.1063/1.3215910
  * Marley et al. (2007a) Marley, M. S., Fortney, J., Seager, S., & Barman, T. 2007a, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 733, doi: 10.48550/arXiv.astro-ph/0602468
  * Marley et al. (2007b) Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007b, ApJ, 655, 541, doi: 10.1086/509759
  * Marley et al. (1995) Marley, M. S., Gómez, P., & Podolak, M. 1995, J. Geophys. Res., 100, 23349, doi: 10.1029/95JE02362
  * Marley & Robinson (2015) Marley, M. S., & Robinson, T. D. 2015, ARA&A, 53, 279, doi: 10.1146/annurev-astro-082214-122522
  * Marley et al. (2012) Marley, M. S., Saumon, D., Cushing, M., et al. 2012, ApJ, 754, 135, doi: 10.1088/0004-637X/754/2/135
  * Marley et al. (2002) Marley, M. S., Seager, S., Saumon, D., et al. 2002, ApJ, 568, 335, doi: 10.1086/338800
  * Marley et al. (2021) Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85, doi: 10.3847/1538-4357/ac141d
  * Marsden et al. (2014) Marsden, S. C., Petit, P., Jeffers, S. V., et al. 2014, MNRAS, 444, 3517, doi: 10.1093/mnras/stu1663
  * Masseron et al. (2016) Masseron, T., Merle, T., & Hawkins, K. 2016, BACCHUS: Brussels Automatic Code for Characterizing High accUracy Spectra, Astrophysics Source Code Library, record ascl:1605.004.  http://ascl.net/1605.004
  * Masseron et al. (2014) Masseron, T., Plez, B., Van Eck, S., et al. 2014, A&A, 571, A47, doi: 10.1051/0004-6361/201423956
  * McQuillan et al. (2014) McQuillan, A., Mazeh, T., & Aigrain, S. 2014, ApJS, 211, 24, doi: 10.1088/0067-0049/211/2/24
  * Mesa et al. (2023) Mesa, D., Gratton, R., Kervella, P., et al. 2023, A&A, 672, A93, doi: 10.1051/0004-6361/202345865
  * Miles et al. (2023) Miles, B. E., Biller, B. A., Patapis, P., et al. 2023, ApJ, 946, L6, doi: 10.3847/2041-8213/acb04a
  * Miller & Fortney (2011) Miller, N., & Fortney, J. J. 2011, ApJ, 736, L29, doi: 10.1088/2041-8205/736/2/L29
  * Mollière et al. (2017) Mollière, P., van Boekel, R., Bouwman, J., et al. 2017, A&A, 600, A10, doi: 10.1051/0004-6361/201629800
  * Mollière et al. (2015) Mollière, P., van Boekel, R., Dullemond, C., Henning, T., & Mordasini, C. 2015, ApJ, 813, 47, doi: 10.1088/0004-637X/813/1/47
  * Mollière et al. (2019) Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67, doi: 10.1051/0004-6361/201935470
  * Mollière et al. (2020) Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131, doi: 10.1051/0004-6361/202038325
  * Mollière et al. (2022) Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74, doi: 10.3847/1538-4357/ac6a56
  * Montalto et al. (2021) Montalto, M., Piotto, G., Marrese, P. M., et al. 2021, A&A, 653, A98, doi: 10.1051/0004-6361/202140717
  * Mordasini et al. (2014) Mordasini, C., Klahr, H., Alibert, Y., Miller, N., & Henning, T. 2014, A&A, 566, A141, doi: 10.1051/0004-6361/201321479
  * Mordasini et al. (2017) Mordasini, C., Marleau, G. D., & Mollière, P. 2017, A&A, 608, A72, doi: 10.1051/0004-6361/201630077
  * Morley et al. (2012) Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012, ApJ, 756, 172, doi: 10.1088/0004-637X/756/2/172
  * Morton (2015) Morton, T. D. 2015, isochrones: Stellar model grid package, Astrophysics Source Code Library, record ascl:1503.010.  http://ascl.net/1503.010
  * Moses et al. (2022) Moses, J. I., Tremblin, P., Venot, O., & Miguel, Y. 2022, Experimental Astronomy, 53, 279, doi: 10.1007/s10686-021-09749-1
  * Mousis et al. (2021) Mousis, O., Lunine, J. I., & Aguichine, A. 2021, ApJ, 918, L23, doi: 10.3847/2041-8213/ac1d50
  * Mukherjee et al. (2023) Mukherjee, S., Batalha, N. E., Fortney, J. J., & Marley, M. S. 2023, ApJ, 942, 71, doi: 10.3847/1538-4357/ac9f48
  * Mukherjee et al. (2022) Mukherjee, S., Fortney, J. J., Batalha, N. E., et al. 2022, ApJ, 938, 107, doi: 10.3847/1538-4357/ac8dfb
  * Nasedkin et al. (2023) Nasedkin, E., Mollière, P., Wang, J., et al. 2023, arXiv e-prints, arXiv:2308.01343, doi: 10.48550/arXiv.2308.01343
  * Nielsen et al. (2017) Nielsen, E. L., Rosa, R. J. D., Rameau, J., et al. 2017, AJ, 154, 218, doi: 10.3847/1538-3881/aa8a69
  * Nielsen et al. (2019) Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13, doi: 10.3847/1538-3881/ab16e9
  * Öberg et al. (2011) Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16, doi: 10.1088/2041-8205/743/1/L16
  * Ohno & Fortney (2022) Ohno, K., & Fortney, J. J. 2022, arXiv e-prints, arXiv:2211.16877, doi: 10.48550/arXiv.2211.16877
  * Ohno & Fortney (2023) —. 2023, ApJ, 946, 18, doi: 10.3847/1538-4357/acafed
  * Oliphant (2006) Oliphant, T. 2006, NumPy: A guide to NumPy, USA: Trelgol Publishing.  http://www.numpy.org/
  * Owen et al. (1999) Owen, T., Mahaffy, P., Niemann, H. B., et al. 1999, Nature, 402, 269, doi: 10.1038/46232
  * Paletou et al. (2015) Paletou, F., Böhm, T., Watson, V., & Trouilhet, J. F. 2015, A&A, 573, A67, doi: 10.1051/0004-6361/201424741
  * Pawellek et al. (2021) Pawellek, N., Wyatt, M., Matrà, L., Kennedy, G., & Yelverton6, B. 2021, MNRAS, 502, 5390, doi: 10.1093/mnras/stab269
  * Pearce et al. (2020) Pearce, L. A., Kraus, A. L., Dupuy, T. J., et al. 2020, ApJ, 894, 115, doi: 10.3847/1538-4357/ab8389
  * Pearce et al. (2022) Pearce, T. D., Launhardt, R., Ostermann, R., et al. 2022, A&A, 659, A135, doi: 10.1051/0004-6361/202142720
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21, doi: 10.1109/MCSE.2007.53
  * Phillips et al. (2020) Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38, doi: 10.1051/0004-6361/201937381
  * Piette & Madhusudhan (2020) Piette, A. A. A., & Madhusudhan, N. 2020, MNRAS, 497, 5136, doi: 10.1093/mnras/staa2289
  * Piso et al. (2016) Piso, A.-M. A., Pegues, J., & Öberg, K. I. 2016, ApJ, 833, 203, doi: 10.3847/1538-4357/833/2/203
  * Plez (2012) Plez, B. 2012, Turbospectrum: Code for spectral synthesis, Astrophysics Source Code Library, record ascl:1205.004.  http://ascl.net/1205.004
  * Polyansky et al. (2018) Polyansky, O. L., Kyuberis, A. A., Zobov, N. F., et al. 2018, MNRAS, 480, 2597, doi: 10.1093/mnras/sty1877
  * Pych (2004) Pych, W. 2004, PASP, 116, 148, doi: 10.1086/381786
  * Radigan et al. (2014) Radigan, J., Lafrenière, D., Jayawardhana, R., & Artigau, E. 2014, ApJ, 793, 75, doi: 10.1088/0004-637X/793/2/75
  * Rajan et al. (2017) Rajan, A., Rameau, J., De Rosa, R. J., et al. 2017, AJ, 154, 10, doi: 10.3847/1538-3881/aa74db
  * Rice & Brewer (2020) Rice, M., & Brewer, J. M. 2020, ApJ, 898, 119, doi: 10.3847/1538-4357/ab9f96
  * Richard et al. (2012) Richard, C., Gordon, I. E., Rothman, L. S., et al. 2012, J. Quant. Spec. Radiat. Transf., 113, 1276, doi: 10.1016/j.jqsrt.2011.11.004
  * Rothman et al. (2010) Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, J. Quant. Spec. Radiat. Transf., 111, 2139, doi: 10.1016/j.jqsrt.2010.05.001
  * Ruffio et al. (2021) Ruffio, J.-B., Konopacky, Q. M., Barman, T., et al. 2021, AJ, 162, 290, doi: 10.3847/1538-3881/ac273a
  * Sahlmann et al. (2011) Sahlmann, J., Ségransan, D., Queloz, D., & Udry, S. 2011, in The Astrophysics of Planetary Systems: Formation, Structure, and Dynamical Evolution, ed. A. Sozzetti, M. G. Lattanzi, & A. P. Boss, Vol. 276, 117–120, doi: 10.1017/S1743921311020047
  * Samland et al. (2017) Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57, doi: 10.1051/0004-6361/201629767
  * Santos et al. (2017) Santos, N. C., Adibekyan, V., Figueira, P., et al. 2017, A&A, 603, A30, doi: 10.1051/0004-6361/201730761
  * Saumon & Marley (2008) Saumon, D., & Marley, M. S. 2008, ApJ, 689, 1327, doi: 10.1086/592734
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schlaufman (2018) Schlaufman, K. C. 2018, ApJ, 853, 37, doi: 10.3847/1538-4357/aa961c
  * Schneider et al. (2014) Schneider, A. C., Cushing, M. C., Kirkpatrick, J. D., et al. 2014, AJ, 147, 34, doi: 10.1088/0004-6256/147/2/34
  * Schneider & Bitsch (2021a) Schneider, A. D., & Bitsch, B. 2021a, A&A, 654, A71, doi: 10.1051/0004-6361/202039640
  * Schneider & Bitsch (2021b) —. 2021b, A&A, 654, A72, doi: 10.1051/0004-6361/202141096
  * Schröder et al. (2009) Schröder, C., Reiners, A., & Schmitt, J. H. M. M. 2009, A&A, 493, 1099, doi: 10.1051/0004-6361:200810377
  * Siess et al. (2000) Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593, doi: 10.48550/arXiv.astro-ph/0003477
  * Skemer et al. (2016) Skemer, A. J., Morley, C. V., Zimmerman, N. T., et al. 2016, ApJ, 817, 166, doi: 10.3847/0004-637X/817/2/166
  * Sousa-Silva et al. (2015) Sousa-Silva, C., Al-Refaie, A. F., Tennyson, J., & Yurchenko, S. N. 2015, MNRAS, 446, 2337, doi: 10.1093/mnras/stu2246
  * Spiegel & Burrows (2012) Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174, doi: 10.1088/0004-637X/745/2/174
  * Suárez & Metchev (2022) Suárez, G., & Metchev, S. 2022, MNRAS, 513, 5701, doi: 10.1093/mnras/stac1205
  * Tagliaferri et al. (1994) Tagliaferri, G., Cutispoto, G., Pallavicini, R., Randich, S., & Pasquini, L. 1994, A&A, 285, 272
  * Takeda et al. (2007) Takeda, G., Ford, E. B., Sills, A., et al. 2007, ApJS, 168, 297, doi: 10.1086/509763
  * Tayar et al. (2022) Tayar, J., Claytor, Z. R., Huber, D., & van Saders, J. 2022, ApJ, 927, 31, doi: 10.3847/1538-4357/ac4bbc
  * Teske et al. (2019) Teske, J. K., Thorngren, D., Fortney, J. J., Hinkel, N., & Brewer, J. M. 2019, AJ, 158, 239, doi: 10.3847/1538-3881/ab4f79
  * Thorngren & Fortney (2019) Thorngren, D., & Fortney, J. J. 2019, ApJ, 874, L31, doi: 10.3847/2041-8213/ab1137
  * Thorngren et al. (2016) Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64, doi: 10.3847/0004-637X/831/1/64
  * Tremblin et al. (2016) Tremblin, P., Amundsen, D. S., Chabrier, G., et al. 2016, ApJ, 817, L19, doi: 10.3847/2041-8205/817/2/L19
  * Tremblin et al. (2019) Tremblin, P., Padioleau, T., Phillips, M. W., et al. 2019, ApJ, 876, 144, doi: 10.3847/1538-4357/ab05db
  * Tsuji et al. (1996) Tsuji, T., Ohnaka, K., Aoki, W., & Nakajima, T. 1996, A&A, 308, L29
  * Valenti & Fischer (2005) Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141, doi: 10.1086/430500
  * Vos et al. (2022) Vos, J. M., Faherty, J. K., Gagné, J., et al. 2022, ApJ, 924, 68, doi: 10.3847/1538-4357/ac4502
  * Vos et al. (2019) Vos, J. M., Biller, B. A., Bonavita, M., et al. 2019, MNRAS, 483, 480, doi: 10.1093/mnras/sty3123
  * Vousden et al. (2016) Vousden, W. D., Farr, W. M., & Mandel, I. 2016, MNRAS, 455, 1919, doi: 10.1093/mnras/stv2422
  * Wahl et al. (2017) Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649, doi: 10.1002/2017GL073160
  * Wang et al. (2020) Wang, J., Wang, J. J., Ma, B., et al. 2020, AJ, 160, 150, doi: 10.3847/1538-3881/ababa7
  * Wang et al. (2022) Wang, J., Kolecki, J. R., Ruffio, J.-B., et al. 2022, AJ, 163, 189, doi: 10.3847/1538-3881/ac56e2
  * Wang et al. (2023) Wang, J., Wang, J. J., Ruffio, J.-B., et al. 2023, AJ, 165, 4, doi: 10.3847/1538-3881/ac9f19
  * Wang et al. (2015) Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015, pyKLIP: PSF Subtraction for Exoplanets and Disks, Astrophysics Source Code Library, record ascl:1506.001.  http://ascl.net/1506.001
  * Wende et al. (2010) Wende, S., Reiners, A., Seifahrt, A., & Bernath, P. F. 2010, A&A, 523, A58, doi: 10.1051/0004-6361/201015220
  * White et al. (2007) White, R. J., Gabor, J. M., & Hillenbrand, L. A. 2007, AJ, 133, 2524, doi: 10.1086/514336
  * Whiteford et al. (2023) Whiteford, N., Glasse, A., Chubb, K. L., et al. 2023, MNRAS, doi: 10.1093/mnras/stad670
  * Williams & Cieza (2011) Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67, doi: 10.1146/annurev-astro-081710-102548
  * Wong et al. (2004) Wong, M. H., Mahaffy, P. R., Atreya, S. K., Niemann, H. B., & Owen, T. C. 2004, Icarus, 171, 153, doi: 10.1016/j.icarus.2004.04.010
  * Wright et al. (2011) Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48, doi: 10.1088/0004-637X/743/1/48
  * Xuan et al. (2022) Xuan, J. W., Wang, J., Ruffio, J.-B., et al. 2022, ApJ, 937, 54, doi: 10.3847/1538-4357/ac8673
  * Yurchenko et al. (2020) Yurchenko, S. N., Mellor, T. M., Freedman, R. S., & Tennyson, J. 2020, MNRAS, 496, 5282, doi: 10.1093/mnras/staa1874
  * Zahnle & Marley (2014) Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41, doi: 10.1088/0004-637X/797/1/41
  * Zalesky et al. (2019) Zalesky, J. A., Line, M. R., Schneider, A. C., & Patience, J. 2019, ApJ, 877, 24, doi: 10.3847/1538-4357/ab16db
  * Zalesky et al. (2022) Zalesky, J. A., Saboi, K., Line, M. R., et al. 2022, ApJ, 936, 44, doi: 10.3847/1538-4357/ac786c
  * Zhang et al. (2021a) Zhang, Y., Snellen, I. A. G., & Mollière, P. 2021a, A&A, 656, A76, doi: 10.1051/0004-6361/202141502
  * Zhang et al. (2021b) Zhang, Y., Snellen, I. A. G., Bohn, A. J., et al. 2021b, Nature, 595, 370, doi: 10.1038/s41586-021-03616-x
  * Zhang et al. (2021c) Zhang, Z., Liu, M. C., Best, W. M. J., Dupuy, T. J., & Siverd, R. J. 2021c, ApJ, 911, 7, doi: 10.3847/1538-4357/abe3fa
  * Zhang et al. (2021d) Zhang, Z., Liu, M. C., Claytor, Z. R., et al. 2021d, ApJ, 916, L11, doi: 10.3847/2041-8213/ac1123
  * Zhang et al. (2021e) Zhang, Z., Liu, M. C., Marley, M. S., Line, M. R., & Best, W. M. J. 2021e, ApJ, 916, 53, doi: 10.3847/1538-4357/abf8b2
  * Zhang et al. (2021f) —. 2021f, ApJ, 921, 95, doi: 10.3847/1538-4357/ac0af7
  * Zhang et al. (2020) Zhang, Z., Liu, M. C., Hermes, J. J., et al. 2020, ApJ, 891, 171, doi: 10.3847/1538-4357/ab765c
  * Zhang et al. (2023) Zhang, Z., Bowler, B. P., Dupuy, T. J., et al. 2023, AJ, 165, 73, doi: 10.3847/1538-3881/aca88c
  * Zhou & Lin (2007) Zhou, J.-L., & Lin, D. N. C. 2007, ApJ, 666, 447, doi: 10.1086/520043
  * Zhou et al. (2016) Zhou, Y., Apai, D., Schneider, G. H., Marley, M. S., & Showman, A. P. 2016, ApJ, 818, 176, doi: 10.3847/0004-637X/818/2/176
  * Zhou et al. (2022) Zhou, Y., Bowler, B. P., Apai, D., et al. 2022, AJ, 164, 239, doi: 10.3847/1538-3881/ac9905
  * Zinn (2021) Zinn, J. C. 2021, AJ, 161, 214, doi: 10.3847/1538-3881/abe936
  * Zúñiga-Fernández et al. (2021) Zúñiga-Fernández, S., Bayo, A., Elliott, P., et al. 2021, A&A, 645, A30, doi: 10.1051/0004-6361/202037830

Figure 16: Our measured stellar parameters of AF Lep A (verticle red line with
shadow) compared to those in the literature (black; sorted by $M_{A}$, and
then [Fe/H]A if $M_{A}$ is not reported). References for labels shown on the y
axis are listed as follows. Tagl94: Tagliaferri et al. (1994); Alle99: Allende
Prieto & Lambert (1999); Felt01: Feltzing et al. (2001); Vale05: Valenti &
Fischer (2005); Ammo06: Ammons et al. (2006); Take07: Takeda et al. (2007);
Casa11: Casagrande et al. (2011); Wrig11: Wright et al. (2011); Pale15:
Paletou et al. (2015); Chan16: Chandler et al. (2016); Luck18: Luck (2018);
Ande19: Anders et al. (2019); Kerv19: Kervella et al. (2019); Rice20: Rice &
Brewer (2020); Mont21: Montalto et al. (2021); Ande22: Anders et al. (2022);
Gaia22: Gaia Collaboration et al. (2022); Pear22: Pearce et al. (2022).
Figure 17: Parameter posteriors for two modes of our derived orbital solution,
corresponding to positive (red) and negative (blue) $\Delta$RVb-A values at
epoch J2024.0. The axes of $\Omega$, $\omega_{\star}$, and $\lambda_{\rm
ref,\star}$ extend beyond the nominal range of [$0^{\circ},360^{\circ}$) so
that each mode of the orbital solution does not split into multiple peaks due
to the modulus of $360^{\circ}$.  Figure 18: Fitted orbits compared to the
observed relative astrometry (top), absolute astrometry (middle), and the host
star’s relative RVs (bottom left). In each panel, we present the observed data
(top) and residuals (bottom) as filled circles, and the best-fit orbital
solution as black solid lines. Predictions of 1000 random orbits are color-
coded by the dynamical mass of AF Lep b, estimated to be $2.7^{+0.6}_{-0.5}$
MJup by our analysis. Figure 19: Results of the retrieval analysis on
$K1/K2/L^{\prime}$ photometry and the Mesa et al. (2023) spectrum of AF Lep b,
with default parameter priors listed in Table 5. The format is the same as
Figure 7. Figure 19: Continued. Results of the retrieval analysis on
$K1/K2/L^{\prime}$ photometry and the Mesa et al. (2023) spectrum of AF Lep b,
with a constrained prior on $M$ and default priors (Table 5) for the remaining
parameters. The format is the same as Figure 7. Figure 19: Continued. Results
of the retrieval analysis on $K1/K2/L^{\prime}$ photometry and the De Rosa et
al. (2023) spectrum of AF Lep b, with default parameter priors listed in Table
5. The format is the same as Figure 7. Figure 19: Continued. Results of the
retrieval analysis on $K1/K2/L^{\prime}$ photometry and the Mesa et al. (2023)
spectrum of AF Lep b, with a constrained prior on $M$ and default priors
(Table 5) for the remaining parameters. The format is the same as Figure 7.
Figure 19: Continued. Results of the retrieval analysis on $K1/K2/L^{\prime}$
photometry and both spectra from Mesa et al. (2023) and De Rosa et al. (2023),
with default parameter priors listed in Table 5. The format is the same as
Figure 7. Figure 19: Continued. Results of the retrieval analysis on
$K1/K2/L^{\prime}$ photometry and both spectra from De Rosa et al. (2023) and
De Rosa et al. (2023), with a constrained prior on $M$ and default priors
(Table 5) for the remaining parameters. The format is the same as Figure 7.
Figure 20: Results of the retrieval analysis on the Mesa et al. (2023)
spectrum of AF Lep b, with default parameter priors listed in Table 5. The
format is the same as Figure 7. Figure 20: Continued. Results of the retrieval
analysis on the Mesa et al. (2023) spectrum of AF Lep b, with a constrained
prior on $M$ and default priors for the remaining parameters (Table 5). The
format is the same as Figure 7. Figure 20: Continued. Results of the retrieval
analysis on the Mesa et al. (2023) spectrum of AF Lep b, with constrained
priors on both $M$ and $R$ and default priors for the remaining parameters
(Table 5). The format is the same as Figure 7. Figure 20: Continued. Results
of the retrieval analysis on the De Rosa et al. (2023) spectrum of AF Lep b,
with default parameter priors listed in Table 5. The format is the same as
Figure 7. Figure 20: Continued. Results of the retrieval analysis on the De
Rosa et al. (2023) spectrum of AF Lep b, with a constrained prior on both $M$
and default priors for the remaining parameters (Table 5). The format is the
same as Figure 7. Figure 20: Continued. Results of the retrieval analysis on
the De Rosa et al. (2023) spectrum of AF Lep b, with constrained priors on
both $M$ and $R$ and default priors for the remaining parameters (Table 5).
The format is the same as Figure 7. Figure 20: Continued. Results of the
retrieval analysis on the both spectra from Mesa et al. (2023) and De Rosa et
al. (2023), with default parameter priors listed in Table 5. The format is the
same as Figure 7. Figure 20: Continued. Results of the retrieval analysis on
$K1/K2/L^{\prime}$ photometry and both spectra from Mesa et al. (2023) and De
Rosa et al. (2023), with a constrained prior on $M$ and default priors for the
remaining parameters (Table 5). The format is the same as Figure 7. Figure 20:
Continued. Results of the retrieval analysis on $K1/K2/L^{\prime}$ photometry
and both spectra from De Rosa et al. (2023) and De Rosa et al. (2023), with
constrained priors on both $M$ and $R$ and default priors for the remaining
parameters (Table 5). The format is the same as Figure 7.

Table 7: Properties of AF Lep b based on evolution models (assuming planet is 10 Myr younger than the host star) Evolution Model |  $T_{\rm eff}$ |  $\log{(L_{\rm bol}/L_{\odot})}$ |  $\log{g}$ |  $R$
---|---|---|---|---
| (K) | (dex) | (dex) | ($R_{\rm Jup}$)
Hot-Start Models
Saumon & Marley (2008): cloudy ($f_{\rm sed}=2$) and [Fe/H]$=0$ | $754^{+72}_{-55}$ | $-5.22^{+0.17}_{-0.15}$ | $3.56^{+0.08}_{-0.07}$ | $1.400^{+0.025}_{-0.022}$
Saumon & Marley (2008): hybrid and [Fe/H]$=0$ | $781^{+106}_{-82}$ | $-5.20^{+0.24}_{-0.21}$ | $3.60^{+0.08}_{-0.07}$ | $1.328^{+0.024}_{-0.018}$
Marley et al. (2021): cloudless and [Fe/H]$=-0.5$ | $744^{+99}_{-82}$ | $-5.31^{+0.23}_{-0.21}$ | $3.62^{+0.08}_{-0.09}$ | $1.288^{+0.018}_{-0.014}$
Marley et al. (2021): cloudless and [Fe/H]$=0$ | $756^{+99}_{-87}$ | $-5.27^{+0.23}_{-0.22}$ | $3.60^{+0.08}_{-0.09}$ | $1.312^{+0.022}_{-0.016}$
Marley et al. (2021): cloudless and [Fe/H]$=+0.5$ | $767^{+95}_{-85}$ | $-5.23^{+0.22}_{-0.22}$ | $3.58^{+0.08}_{-0.09}$ | $1.339^{+0.026}_{-0.020}$
Spiegel & Burrows (2012): cloudless and [Fe/H]$=0$ | $742^{+101}_{-89}$ | $-5.29^{+0.23}_{-0.23}$ | $3.59^{+0.09}_{-0.09}$ | $1.329^{+0.023}_{-0.017}$
Spiegel & Burrows (2012): hybrid and [Fe/H]$=+0.5$ | $723^{+98}_{-86}$ | $-5.34^{+0.23}_{-0.23}$ | $3.59^{+0.09}_{-0.09}$ | $1.329^{+0.023}_{-0.017}$
Phillips et al. (2020): cloudless and [Fe/H]$=0$ | $760^{+99}_{-91}$ | $-5.26^{+0.23}_{-0.23}$ | $3.60^{+0.08}_{-0.09}$ | $1.318^{+0.027}_{-0.020}$
Baraffe et al. (2003): cloudless and [Fe/H]$=0$ | $791^{+101}_{-90}$ | $-5.20^{+0.22}_{-0.22}$ | $3.60^{+0.08}_{-0.09}$ | $1.317^{+0.023}_{-0.020}$
Cold-Start Models
Marley et al. (2021): cloudless and [Fe/H]$=-0.5$ | $665^{+41}_{-44}$ | $-5.52^{+0.11}_{-0.12}$ | $3.62^{+0.09}_{-0.10}$ | $1.278^{+0.013}_{-0.010}$
Spiegel & Burrows (2012): cloudless and [Fe/H]$=0$ | $520^{+21}_{-22}$ | $-5.96^{+0.07}_{-0.08}$ | $3.64^{+0.10}_{-0.10}$ | $1.253^{+0.012}_{-0.010}$
Spiegel & Burrows (2012): hybrid and [Fe/H]$=+0.5$ | $528^{+24}_{-26}$ | $-5.93^{+0.08}_{-0.09}$ | $3.64^{+0.10}_{-0.10}$ | $1.253^{+0.012}_{-0.010}$


# MPC of Uncertain Nonlinear Systems with Meta-Learning for Fast Adaptation of
Neural Predictive Models

Jiaqi Yan, Ankush Chakrabarty, Alisa Rupenyan, and John Lygeros  J. Yan and J.
Lygeros are with Automatic Control Laboratory, ETH Zurich, Switzerland
({jiayan,jlygeros}@ethz.ch). A. Chakrabarty is with Mitsubishi Electric
Research Laboratories, Cambridge, MA 02139, USA (achakrabarty@ieee.org). A.
Rupenyan is with ZHAW Centre for Artificial Intelligence, Zurich University of
Applied Sciences, Switzerland (alisa.rupenyan@zhaw.ch). This work is supported
by the Swiss National Science Foundation supported this work through NCCR
Automation under Grant agreement 51NF40_180545.

###### Abstract

In this paper, we consider the problem of reference tracking in uncertain
nonlinear systems. A neural State-Space Model (NSSM) is used to approximate
the nonlinear system, where a deep encoder network learns the nonlinearity
from data, and a state-space component captures the temporal relationship.
This transforms the nonlinear system into a linear system in a latent space,
enabling the application of model predictive control (MPC) to determine
effective control actions. Our objective is to design the optimal controller
using limited data from the target system (the system of interest). To this
end, we employ an implicit model-agnostic meta-learning (iMAML) framework that
leverages information from source systems (systems that share similarities
with the target system) to expedite training in the target system and enhance
its control performance. The framework consists of two phases: the (offine)
meta-training phase learns a aggregated NSSM using data from source systems,
and the (online) meta-inference phase quickly adapts this aggregated model to
the target system using only a few data points and few online training
iterations, based on local loss function gradients. The iMAML algorithm
exploits the implicit function theorem to exactly compute the gradient during
training, without relying on the entire optimization path. By focusing solely
on the optimal solution, rather than the path, we can meta-train with less
storage complexity and fewer approximations than other contemporary meta-
learning algorithms. We demonstrate through numerical examples that our
proposed method can yield accurate predictive models by adaptation, resulting
in a downstream MPC that outperforms several baselines.

##  I Introduction

Optimal control for unknown nonlinear systems has been a long-standing
challenge in various applications such as robotics, manufacturing systems, and
so on [1, 2, 3]. Different learning-based approaches have been explored, such
as reinforcement learning, adaptive dynamic programming, and stochastic
optimization-based control [4, 5, 6]. Recently, the neural state-space model
(NSSM) has gained popularity as an effective tool for tackling nonlinear
control problems [7, 8, 9]. This model enhances traditional state-space models
by employing neural networks to capture complex and nonlinear relationships
within the system. Specifically, the neural network component empowers the
model to learn from data and transform the original state-space into a latent
space where the dynamics are reasonably approximated by a linear system [10].
This architecture is grounded in Koopman operator theory, demonstrating that
any nonlinear system can be lifted to an infinite-dimensional linear system in
a latent space [11]. For the sake of tractable computations, the latent state
of NSSMs is chosen to be of finite size [12].

Although NSSMs show good performance in approximating nonlinear systems, their
training usually requires a large dataset. In situations where data is scarce,
constructing an accurate model and developing an efficient learning-based
controller becomes challenging. This difficulty often arises when collecting
data from the system of interest, referred to as the target system, is arduous
or expensive. In response to this challenge, recent efforts, while not
specifically within the domain of NSSMs, aim to leverage data collected from
source systems that share similarities with the target system [13, 14, 15]. As
an example, the target system could be a physical system, while the source
system may be a numerical model or a digital twin of the target system, or a
system with modified physical parameters that is easily accessible, allowing
for the collection of large amounts of data. The objective is to use this
shared dataset to pre-train a model, then fine-tune it on the target system
with limited data.

In this paper, we extend the approach of leveraging knowledge acquired from
similar systems to expedite training of the target system in the context of
NSSM. Different from existing works [13, 15] that use a transfer learning
framework, we propose a method by using meta-learning. The distinction between
transfer learning and meta-learning is elucidated in [16], and differences in
empirical performance between the two classes of methods have been highlighted
in [7].

Meta-learning has recently gained widespread application in control and
robotics, addressing challenges such as system identification, parameter
estimation, and adaptive control [17, 7, 18, 19, 20]. In most of these works,
model-agnostic meta-learning (MAML) based solutions are developed, which
involves solving a bilevel optimization problem in the meta-training stage.
The inner loop of MAML involves calculating and propagating derivatives of the
training loss function along the full optimization path (which is the same as
the number of inner-loop updates/steps). The outer loop then accumulates
information from these inner-loop updates and computes an outer-loop gradient
direction. By repeating these two loops, the MAML learner is expected to
asymptotically converge to a set of neural weights from which rapid adaptation
is possible to any task within the distribution of training tasks. Choosing a
large number of inner-loop steps, while useful to understand good outer-loop
directions, incurs high memory complexity, and is not practical for a large
number of inner-loop adaptation steps. Therefore, existing approaches
typically propose approximations to obtain a more efficient solution or to
allow for a large number of inner-loop steps. Unfortunately, these
approximations are sometimes oversimplifications or valid for specific
activation functions (e.g., first-order MAML exploits ReLU properties), and
thus, can harm quality. One particular approximation that utilizes the
implicit function theorem is called implicit MAML (iMAML), that enables
computing the asymptotic inner-loop gradient efficiently [21]. By doing this,
the meta-learner can update the outer-loop based on directions obtained by
solving the inner-loop over a large number of iterations without requiring
high memory, as the method is agnostic to the path taken in the inner-loop and
only uses the final gradient direction. As a result, better adaptation
performance can be expected.

In this paper, we solve the optimal tracking problem in unknown nonlinear
systems by leveraging the iMAML algorithm for system identification, and using
this identified predictive model within a receding horizon framework. This is
summarized herein:

(i) To efficiently utilize limited samples from the target system, we propose
a meta-learning framework that customizes iMAML for adapting neural predictive
models. The framework pre-trains an aggregated model from source systems, and
then fine-tunes it with a small target dataset and few adaptation steps. (ii)
We incorporate the adapted models into MPC, determining optimal inputs for the
target system. (iii) Numerical examples show that our method outperforms
models trained solely by MAML or with only target system data, along with
other baselines.

Organization: Section II introduces the problem of interest. In Section III,
we propose the meta-learning-based MPC to solve the problem by training the
NSSM with data from source systems. The performance of our algorithm is
evaluated through numerical examples in Section IV. Finally, Section V
concludes the paper.

##  II Problem Formulation

We consider a family of parameterized discrete-time nonlinear systems of the
form

$\begin{split}x_{t+1}&=f(x_{t},u_{t},\theta_{f}),\\\
y_{t}&=g(x_{t},u_{t},\theta_{g}),\end{split}$ (1)

where $x\in\mathbb{R}^{n}$ and $y\in\mathbb{R}^{m}$ are respectively system
state and output. Moreover, $u\in\mathcal{C}\subseteq\mathbb{R}^{p}$ is the
input with $\mathcal{C}$ being a compact, convex set. Both $f$ and $g$ are
unknown nonlinear functions, and
$\theta:=[\theta_{f};\theta_{g}]\in\mathbb{R}^{w}$ denotes a vector of unknown
parameters. Despite the unknown dynamics, our objective is to design a
controller to asymptotically track some reference signal $\bar{y}$ in the
sense that

$\lim_{t\to\infty}y_{t}=\bar{y}_{t}.$ (2)

Since we aim to solve the problem by using a meta-learning framework, we
assume that the target system under consideration is similar to some systems,
from which we have archival data.

###  II-A Neural State-Space Model (NSSM)

For the unknown system described by (1), we construct a neural state-space
model (NSSM) that approximates the dynamics of the target system (the system
of interest) parameterized by $\theta=\theta^{*}$. Based on the NSSM, the
optimal controller can be designed using a model-based approach. The true
value of $\theta^{*}$ is unknown, but it is drawn from a distribution

$\theta^{*}\sim\Theta.$ (3)

The distribution $\Theta$ is often informed by domain experience, e.g., a
particular parameter may have some interpretation from physics that defines a
range of admissible values, and one may sample from a uniform distribution
over this range. Therefore, it is not impractical to assume we know $\Theta$,
but it is not necessarily for our proposed method: we instead assume that
$\Theta$ can be reliably sampled from even though its distributional form is
unknown.

Let $U(\theta^{*},T^{*})$ and $Y(\theta^{*},T^{*})$ represent the input and
output trajectories generated from the target system over a horizon $T^{*}$
and construct a target dataset

$\mathcal{D}_{\text{target}}:=\\{U(\theta^{*},T^{*}),Y(\theta^{*},T^{*})\\}.$
(4)

Using this dataset, one can train an NSSM of the form

$\displaystyle z_{t}$ $\displaystyle=f_{\text{enc}}(U_{t-H+1:t},Y_{t-H+1:t}),$
(5a) $\displaystyle z_{t+1}$ $\displaystyle=A_{z}z_{t}+B_{z}u_{t+1},$ (5b)
$\displaystyle\hat{y}_{t}$ $\displaystyle=C_{z}z_{t}.$ (5c)

where $z\in\mathbb{R}^{n_{z}}$ is the latent state learned by the encoder
network $f_{\text{enc}}(\cdot)$ from historical input and output data:

$\begin{split}U_{t-H+1:t}&:=\\{u_{t-H+1},u_{u-H+2},\cdots,u_{t}\\}\subseteq\mathcal{D}_{\text{target}},\\\
Y_{t-H+1:t}&:=\\{y_{t-H+1},y_{t-H+2},\cdots,y_{t}\\}\subseteq\mathcal{D}_{\text{target}}.\end{split}$
(6)

Here, the latent state dimension $n_{z}\in\mathbb{N}$ and window length
$H\in\mathbb{N}$ are user-defined parameters; these are standard
hyperparameters in classical system identification as well. Also,
$\hat{y}_{t}\in\mathbb{R}^{m}$ is the predicted output, which we require to be
close to the true output $y_{t}$ after NSSM training [7].

Note that the NSSM (5) consists of two components: the deep latent encoder
$f_{\text{enc}}$ in (5a) captures the nonlinearity of the system, while the
state-space component in (5b)–(5c) models temporal dependencies in the system.
In practical implementations, $f_{\text{enc}}$ is trained from standard deep
learning. On the other hand, the state transition in (5b) and (5c) can be
modeled by a recurrent neural network (RNN). The RNN layer takes $u$ and $z$
as the current input and hidden state, respectively, and outputs the next
prediction $\hat{y}$.

Training (5) involves optimizing the weights of the encoder network
$f_{\text{enc}}(\cdot)$, and the elements of the linear decoders
$A_{z},B_{z},C_{z}$ and $D_{z}$; let $\omega$ be the set of these parameters.
The prediction performance of (5) is evaluated as follows.

  1. 1.

We first construct a dataset $\mathcal{D}$ from $\mathcal{D}_{\text{target}}$
by collecting trajectories of length $H+T$:

$\begin{split}U_{t-H+1:t+T}&:=\\{u_{t-H+1},u_{u-H+2},\cdots,u_{t+T}\\},\\\
Y_{t-H+1:t+T}&:=\\{y_{t-H+1},y_{t-H+2},\cdots,y_{t+T}\\}.\end{split}$

  2. 2.

With $U_{t-H+1:t}$ and $Y_{t-H+1:t}$, we then obtain the latent state $z_{t}$
using (5a).

  3. 3.

For a prediction horizon of $T$ and with input $U_{t+1:t+T}$, let us
recursively compute

$\hat{Y}_{t+1:t+T}:=\\{\hat{y}_{t+1},\hat{y}_{t+2},\cdots,\hat{y}_{t+T}\\}$

from (5b) and (5c).

  4. 4.

The performance of the NSSM is evaluated through the following loss function

$\ell_{\text{SSM}}=\frac{1}{T}||Y_{t+1:t+T}-\hat{Y}_{t+1:t+T}||^{2}.$

Obviously, $\ell_{\text{SSM}}$ is determined by both the dataset and the
parameters of (5). To make it clearer, we explicitly define
$\ell_{\text{SSM}}(\operatorname*{\mathcal{D}};\omega)$ as the loss function
evaluated on the dataset $\operatorname*{\mathcal{D}}$ and the parameter
$\omega$.

###  II-B NSSM-enabled MPC

Once trained, the NSSM can be used for predictive control tasks, such as
within an MPC framework. To leverage this model to calculate the optimal input
for the reference tracking, we rewrite (5b) and (5c) in a compact form

$s_{t+1}=As_{t}+Bu_{t+1},$ (7)

where

$s_{t+1}:=\begin{bmatrix}z_{t+1}\\\
\hat{y}_{t+1}\end{bmatrix},\;A:=\begin{bmatrix}A_{z}&0\\\
C_{z}A_{z}&0\end{bmatrix},\;B:=\begin{bmatrix}B_{z}\\\
C_{z}B_{z}\end{bmatrix}.$

With this linearized model, we propose to track the reference based on model
predictive control (MPC). Notice that the choice of the controller is not
limiting for the approach. Here we focus on MPC as it is a suitable approach
for trajectory tracking under constraints, as demonstrated in ([22]). It
solves the following trajectory tracking problem:

$\begin{split}\min_{u}\;&\tilde{y}_{t+N\mid t}^{T}\cdot
P\cdot\tilde{y}_{t+N\mid t}\\\ &+\sum_{k=0}^{N-1}(\tilde{y}_{t+k\mid
t}^{T}\cdot Q\cdot\tilde{y}_{t+k\mid t}+\Delta u_{t+k\mid t}^{T}\cdot
R\cdot\Delta u_{t+k\mid t})\\\ \text{ s.t. }\;&\eqref{eqn:SSM}\text{ and
}u_{t+k\mid t}\in\mathcal{C},\;k=0,\ldots,N-1,\end{split}$ (8)

where for $k=0,\ldots,N-1$, we define

$\begin{split}\tilde{y}_{t+k\mid t}&\triangleq\hat{y}_{t+k\mid
t}-\bar{y}_{t+k},\\\ \Delta u_{t+k\mid t}&\triangleq u_{t+k+1\mid
t}-u_{t+k\mid t}.\end{split}$

Moreover, $Q\geq 0$ and $R>0$ are design parameters kept fixed, and $P$ is the
solution to the discrete-time algebraic Ricatti equation

$P=A^{T}\left(P-PB\left(R+B^{T}PB\right)^{-1}B^{T}P\right)A+Q.$

Note that in (5), the nonlinearity is just used to initialize an otherwise
linear prediction model. Hence the MPC optimization problem is convex and
computationally tractable. Also notice that the formulation (8) is not unique.
Other formulations (with different terminal costs, add output constraints,
etc.) can also be used here.

The MPC finds the optimal controller that minimizes the cost function in (8),
thereby directing the predicted output $\hat{y}$ towards the reference signal
$\bar{y}$. When the NSSM accurately captures the nonlinear dynamics in (1), we
achieve the reference tracking as claimed in (2).

###  II-C Meta-learning framework

Collecting data is often expensive in many applications. Therefore,
$\mathcal{D}_{\text{target}}$ is typically of limited size. Training with
limited data can lead to poor prediction performance of NSSMs, which can
adversely affect predictive control performance. To address this challenge, we
propose a meta-learning framework. The framework is designed to acquire an
NSSM, usually offline, that leverages data from multiple source systems that
share similarities with the target system.

For instance, the source system could be a numerical model or a digital twin
of the target system, and is easily accessible. It could also be a system with
modified physical parameters (e.g., different load, friction parameters, etc.)
This aggregated model enables quick online adaptation to the target system
through few-shot learning, requiring only a small number of data points.

As another motivating example, suppose that we successfully control some
manufacturing process governed by a nonlinear system, and have access to data
from this system. However, the manufacturing requirements change, which
corresponds to a change in the physical parameters in the system, and the
controller performance needs to be adjusted. By using meta-learning, we
propose a method that achieves this adjustment using very few actual
measurements from the new system (the target system), by shifting the data
requirements to multiple previously observed systems (the source systems).

Assume that we have access to a dataset containing input and output
trajectories generated by $N_{s}$ different source systems in the form of (1),
where the $k$-th system is parameterized by some vector $\theta_{k}$ sampled
from the distribution $\Theta$. It is important to know that $\theta_{k}$ need
not to be known. We represent the source dataset as

$\mathcal{D}_{\text{source}}=\\{(U(\theta^{k},T^{k}),Y(\theta^{k},T^{k}))\\}_{k=1}^{N_{s}}:=\\{\mathcal{D}^{k}\\}_{k=1}^{N_{s}},$
(9)

where $T^{k}$ is the length of trajectory produced by the $k$-th system, and
$\theta_{k}\sim\Theta$ for each $k\in\\{1,\cdots,N_{s}\\}$. Our objective is
to learn an aggregated NSSM by leveraging the source data
$\mathcal{D}_{\text{source}}$. This aggregated model should have the ability
to be quickly adapted to other systems only using small datasets and without
explicitly estimating $\theta$. Later, we discuss how to update the dataset by
incorporating the trajectories produced by MPC on the source systems.

##  III Meta-Learning-Based MPC Design

We tailor the meta-learning algorithm iMAML [21] to address the problem of
reference tracking on the target system. The framework consists of two phases:
the meta-training phase learns the aggregated NSSM using
$\mathcal{D}_{\text{source}}$, while the meta-inference phase adapts this
model to the target system using $\mathcal{D}_{\text{target}}$.

###  III-A Bi-level optimization problem in meta-learning

The phase of meta-training involves two loops: the outer-loop updates the
parameters of the NSSM based on the performance across multiple source systems
and the inner-loop adapts the model to individual source systems. This enables
the meta-learning algorithm to learn from multiple systems and generalize well
to new datasets.

Figure 1: Information flow in the meta-training phase.

Specifically, as shown in Fig. 1, the outer loop provides an aggregated
representation $\omega$ of the NSSM (5). In each iteration of the outer loop,
we sample a batch of trajectories
$\\{\operatorname*{\mathcal{D}}^{b}\\}_{k=1}^{B}$ from the source dataset
$\mathcal{D}{\text{source}}$. Then the inner-loop trains individual weights
$\omega^{b}$ using the aggregated representation and the source data. After
this, the outer loop updates the aggregated representation $\omega$ using
these individual values. The process is repeated until convergence.

For each $\operatorname*{\mathcal{D}}^{b}$, we partition it into a training
set $\operatorname*{\mathcal{D}}^{b}_{\text{tr}}$ and a testing set
$\operatorname*{\mathcal{D}}^{b}_{\text{test}}$. The training set is utilized
to adapt the aggregated representation $\omega$ and produce the task-specific
parameter $\omega^{b}$, tailored to the system (1) parameterized by
$\theta^{b}$. To be specific, within the context of iMAML, the inner-loop
solves the following optimization problem:

$\omega^{b}(\omega):=\operatorname*{argmin}_{\psi}\;\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{tr}};\psi)+\frac{\gamma}{2}||\psi-\omega||^{2}.$
(10)

Here, $\omega$ is the weight of the aggregated NSSM, and the regularization
term $||\psi-\omega||^{2}$ encourages the finding of an optimal solution
within a (ideally, small) neighborhood of $\omega$, and $\gamma>0$ controls
the regularization strength. As we will shown later in Sections III-B and
III-D, this regularization enables us to express the outer-loop gradient term
in closed-form by exploiting the implicit function theorem. By doing so, the
gradient can be calculated exactly, without relying on the inner-loop
optimization path but solely on the optimal solution. As a result, there is no
need to store the full optimization path or approximate the gradient, as
required in MAML [23]. This approach reduces memory consumption and enhances
predictive performance, compared to MAML whose memory complexity is high for a
large number of inner-loop steps.

The testing set $\mathcal{D}^{b}_{\text{test}}$ is employed in the outer loop
to update the aggregated representation $\omega$. This is achieved by
evaluating the prediction performance of the NSSM across multiple systems, as
defined below

$\min_{\omega}\;\frac{1}{B}\sum_{b=1}^{B}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{test}};\omega^{b}(\omega)).$
(11)

As shown in (10), $\omega^{b}(\omega)$ is the optimal task-specific parameter
adapted from $\omega$. Therefore, the objective of the outer loop is to learn
a set of parameters $\omega$ that can produce good task-specific parameters
after adaptation. For simplicity, we will denote $\omega^{b}(\omega)$ as
$\omega^{b}$ in the rest of this paper.

###  III-B Approximation solutions to the bi-level optimization problem

To efficiently solve the bi-level optimization problem (10) and (11), the
solution of (11) can be approached using the gradient decent algorithm:

$\omega\leftarrow\omega-\beta_{\text{out
}}\frac{1}{B}\sum_{b=1}^{B}\nabla_{\omega}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{test}};\omega^{b}),$
(12)

where $\beta_{\text{out }}$ denotes the learning rate for the outer loop.
Considering (10), $\omega^{b}$ also depends on $\omega$. Therefore, by
applying the chain rule, we expand the gradient term in (12) as

$\nabla_{\omega}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{test}};\omega^{b})=\frac{d\omega^{b}}{d\omega}\cdot
P^{b}.$ (13)

Here,

$P^{b}:=\nabla_{\psi}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{test}};\psi)|_{\psi=\omega^{b}},$

which is the value of
$\nabla_{\psi}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{test}};\psi)$
evaluated at $\psi=\omega^{b}$. Since $P^{b}$ can be easily obtained via back
propagation on the NSSM, we focus on calculating
$\frac{d\omega^{b}}{d\omega}$.

Recall that $\omega^{b}$ is obtained by solving the inner optimization problem
(10), typically through iterative algorithms like gradient descent. Thus, one
approach to compute $\frac{d\omega^{b}}{d\omega}$ is propagating derivatives
throughout the iterative process in the inner-loop. As this requires the full
path of optimization to be stored in the memory, the approach becomes
intractable if the number of iterative steps is large, or if the inner-loop
procedure is non-differentiable. To address this challenge, the following
lemma introduces an alternate method to calculate this gradient, which does
not rely on the optimization path, but only on the optimal solution
$\omega^{b}$.

######  Lemma 1 ([21])

Define

$Q^{b}:=I+\frac{1}{\gamma}\nabla_{\psi}^{2}\ell_{\text{SSM}}(\mathcal{D}^{b}_{\text{tr}};\psi)|_{\psi=\omega^{b}},$
(14)

where
$\nabla_{\psi}^{2}\ell_{\text{SSM}}(\mathcal{D}^{b}_{\text{tr}};\psi)|_{\psi=\omega^{b}}$
is value of the Hessian matrix
$\nabla_{\psi}^{2}\ell_{\text{SSM}}(\mathcal{D}^{b}_{\text{tr}};\psi)$
evaluated at $\psi=\omega^{b}$. If $Q^{b}$ is invertible, then

$\frac{d\omega^{b}}{d\omega}=(Q^{b})^{-1}.$ (15)

Combining (13) and (15), we obtain

$\nabla_{\omega}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{test}};\omega^{b})=(Q^{b})^{-1}P^{b}.$
(16)

As $Q^{b}$ solely depends on $\omega^{b}$, Lemma 1 offers a solution for
computing $\frac{d\omega^{b}}{d\omega}$ with reduced memory requirements.
Nevertheless, obtaining $\omega^{b}$, the exact solution to the inner problem,
involves solving (10) until convergence, a task often impractical in practice.
Therefore, in this paper, we resort to employing the gradient descent
algorithm with a finite number of steps to obtain an approximate solution.

On the other hand, calculating the inverse of $Q^{b}$ can also be expensive,
especially for large NSSMs. In order to tackle this issue, we notice that
$(Q^{b})^{-1}P^{b}$ is the solution to the following problem:

$\min_{\phi}\;\phi^{T}Q^{b}\phi-\phi^{T}P^{b}.$ (17)

Therefore, in practice, we can also approximate $(Q^{b})^{-1}P^{b}$ by solving
(17) with a finite number of iterations, e.g., update the solution by using a
few gradient steps. By doing so, one avoids the need to explicitly compute the
matrix inverse.

Finally, in order to ensure the compatibility with MPC, we update the source
dataset online. After adaptation in each inner loop, we determine the optimal
input $u^{b}_{*}$ by solving the MPC problem (8) using the NSSM parameterized
by $\omega^{b}$. To balance exploration and exploitation, we let

$u^{b}=u^{b}_{*}+u^{b}_{0},$ (18)

where $u^{b}_{0}$ is a random variable drawn from a Gaussian distribution,
i.e., $u^{b}_{0}\sim\mathcal{N}(0,\Sigma^{b})$. Injecting $u^{b}$ into the
system (1) (parameterized by $\theta^{b}$) yields an output $y^{b}$. We then
update $\mathcal{D}_{\text{source}}$ by incorporating the most recent
trajectory $(u^{b},y^{b})$. The complete meta-training algorithm is summarized
in Algorithm 1.

Algorithm 1 Meta-training phase

0: $\omega\leftarrow$ randomly initialize the parameters of the NSSM

0: $\mathcal{D}_{\text{source }}\leftarrow$ initial source dataset

0: Regularization strength $\gamma$ and learning rate $\beta_{\text{out}}$

while not converge do $\triangleright$ outer-loop

Sample batch $\\{\operatorname*{\mathcal{D}}^{b}\\}_{b=1}^{B}$ from
$\mathcal{D}_{\text{source }}$

for $b=1:B$ do $\triangleright$ inner-loop

Partition $\operatorname*{\mathcal{D}}^{b}$ into $\mathcal{D}^{b}_{\text{tr}}$
and $\mathcal{D}^{b}_{\text{test}}$

$\omega^{b}\leftarrow$ solving (10) by gradient descent

$g^{b}\leftarrow$ solving (17) by gradient descent

$u^{b}\leftarrow$ the new input computed by (8) and (18) with the NSSM
parameterized by $\omega^{b}$

$y^{b}\leftarrow$ simulating the system (1) using $u^{b}$ with
$\theta=\theta^{b}$

Update $\mathcal{D}_{\text{source}}$ by adding $(u^{b},y^{b})$

end for

$\omega\leftarrow\omega-\beta_{\text{out }}\frac{1}{B}\sum_{b=1}^{B}g^{b}$

end while

$\omega_{\infty}\leftarrow$ parameters obtained after meta-training phase

###  III-C Meta-inference

After the meta-learning phase, we acquire the aggregated NSSM parameterized by
$\omega_{\infty}$. The subsequent meta-inference phase adapts this aggregated
model to accurately capture the characteristics of the target system and
achieve the optimal tracking.

Algorithm 2 Meta-inference phase with tracking MPC

0: $\omega_{\infty}\leftarrow$ parameters obtained after meta-training phase

0: $\mathcal{D}_{\text{target }}\leftarrow$ initial target dataset

0: $\gamma$

$\omega^{*}\leftarrow$ solving (19) by gradient descent $\triangleright$
adaptation

while not done do $\triangleright$ reference tracking

$u^{*}\leftarrow$ the optimal solution of (8) based on the NSSM parameterized
by $\omega^{*}$

Inject $u^{*}$ to the target system

end while

Algorithm 2 outlines the meta-inference process; note that no outer-loop is
involved in this phase. We only perform an adaptation step similar to (10) to
derive the specific NSSM for the target system by using the limited dataset
$\mathcal{D}_{\text{target}}$. This is expressed as

$\omega^{*}:=\operatorname*{argmin}_{\psi}\;\ell_{\mathrm{SSM}}(\mathcal{D}_{\text{target}};\psi)+\frac{\gamma}{2}||\psi-\omega_{\infty}||^{2}.$
(19)

In practice, we approximate the optimal solution using only a few gradient
steps. This corresponds to very few “trials” to correctly adapt the
parameters.

Now, with this NSSM tailored for the target system, i.e., the NSSM
parameterized with $\omega^{*}$, we can achieve the reference tracking by
finding the optimal input through MPC (8).

###  III-D Comparison with MAML

Algorithm 1 is developed based on the iMAML algorithm [21]. Alternatively, one
could also consider here the MAML algorithm [23]. Unlike (10), the inner loop
of the MAML-based NSSM approximates the solution of

$\omega^{b}:=\operatorname*{argmin}_{\psi}\;\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{tr}};\psi).$
(20)

Namely, setting $\gamma$ in (10) as $0$. In practice, the algorithm performs
the adaptation through the gradient descent:

$\begin{split}\omega_{0}^{b}&=\omega,\\\
\omega_{m}^{b}&=\omega_{m-1}^{b}-\beta_{\text{in
}}\nabla_{\omega_{m-1}^{b}}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{tr}};\omega_{m-1}^{b}),\;i=1,\cdots,M,\\\
\omega^{b}&=\omega_{M}^{b},\end{split}$ (21)

where $M$ is the number of inner-loop iterations and $\beta_{\text{in }}$ is
the inner-loop learning rate. The outer-loop of MAML aims to solve the same
problem as (11). However, since a different cost function is optimized in
(20), $\frac{d\omega^{b}}{d\omega}$ cannot be directly calculated as in (15).
The exact calculation involves propagating derivatives along the optimization
path, which is not memory-efficient in practice. To tackle this, existing
approaches typically resort to approximating (12) in the outer loop update as:

$\omega\leftarrow\omega-\beta_{\text{out
}}\frac{1}{B}\sum_{b=1}^{B}\nabla_{\omega^{b}}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{test}};\omega^{b}),$
(22)

where it is assumed that

$\nabla_{\omega}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{test}};\omega^{b})\approx\nabla_{\omega^{b}}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{test}};\omega^{b}).$
(23)

This approximation results in a tractable solution since
$\nabla_{\omega^{b}}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{test}};\omega^{b})$
can be readily obtained. However, (23) may harm prediction performance of the
NSSMs. In contrast, by leveraging the iMAML, our proposed algorithm can
calculate the gradient without relying on the optimization path, see (16).
Therefore, there is no need to store the full optimization path or approximate
the gradient as required in MAML. Consequently, the proposed approach reduces
memory consumption and enhances control performance. In Section IV, we will
compare the performance of our algorithm and the MAML-based one through
numerical examples.

##  IV Simulation

In this section, we illustrate the performance of our meta-learned MPC through
some numerical examples.

###  IV-A Van der Pol oscillators

We first consider a family of van der Pol oscillators. The dynamics of each
oscillator is given by

$\begin{split}\dot{x}_{1}&=x_{2},\;\dot{x}_{2}=\theta
x_{2}\left(1-x_{1}^{2}\right)-x_{1}+u,\\\ y&=x,\end{split}$ (24)

where $\theta$ is the unknown damping ratio. By using meta-learning, we want
the approach to work for various different damping ratios and to adapt given
the data of several of them. The source and target systems are generated by
sampling $\theta$ from a Gaussian distribution $\theta\sim\mathcal{N}(0,1)$.

We select the dimension of the latent state $z$ as $n_{z}=5$. The encoder
$f_{\text{enc}}$ is formed by a neural network comprising one input layer, one
hidden layer, and one output layer, with each layer containing $128$ neurons
and activated by rectified linear units (ReLUs). Notice that we can afford to
use a shallow network here due to the explicit incorporation of the adaptation
mechanism into the meta-training process. However, if we were required to
learn a aggregated model for the entire distribution $\Theta$, we would need
more depth and a larger-dimensional latent vector. The matrices in the SSM
$A_{z},B_{z}$ and $C_{z}$, are randomly initialized. The experiments are
repeated for $10$ different random seeds.

First, we evaluate the prediction performance of the NSSM during the meta-
training phase. We compare the iMAML-based Algorithm 1, with the MAML-based
approach as outlined in Section III-D.

The comparison is depicted in Fig. 2, where the $x$-axis refers to the number
of outer loop updates, that is, (12) and (22) respectively in the two
algorithms. In each outer loop, we sample a batch of $B=16$ datasets from
$\mathcal{D}_{\text{source}}$. Moreover, the performance is evaluated by
taking the logarithm of
$1/B\sum_{b=1}^{B}\ell_{\mathrm{SSM}}(\mathcal{D}^{b}_{\text{test}};\omega^{b})$,
showing the prediction performance of the NSSMs. It is observed that Algorithm
1 outperforms the MAML-based approach throughout the training stage since it
computes the derivative $\frac{d\omega^{b}}{d\omega}$ more accurately and thus
requires less iterations for the same loss, as discussed in Section III-D.
Therefore, it results in a smaller approximation error of the NSSM.

Figure 2: Comparison of the prediction performance in the meta-training phase.
The lines and shaded regions respectively represent the mean and standard
deviation across $10$ runs.

To evaluate the tracking performance of the proposed algorithm on the target
system, we compare against a few approaches, all of which have the same NSSM
architecture but are initialized differently:

  1. 1.

iMAML: The parameters of the NSSM are set with those obtained from the meta-
training phase of the iMAML-based algorithm, that is, $\omega_{\infty}$ output
by Algorithm 1.

  2. 2.

MAML: The parameters of the NSSM are set with those obtained from the meta-
training phase of the MAML-based algorithm as outlined in Section III-D. While
Fig. 2 shows that iMAML yields less prediction errors over MAML during the
meta-training phase, we would like to investigate whether this leads to a
better control performance on the target system.

  3. 3.

Supervised learning: Since there is no meta-training phase, the parameters are
chosen randomly. This comparison with Supervised learning highlights whether
meta-learning can improve our performance by leveraging the knowledge from
source systems.

Starting from these different initial parameters, we train the NSSMs using
data from the target system $\operatorname*{\mathcal{D}}_{\text{target}}$,
which only contains $300$ data points. Specifically, for iMAML, this refers to
run the adaptation step in Algorithm 2. On the other hand, iMAML and
supervised learning perform the following update

$\begin{split}\omega_{m}^{*}&=\omega_{m-1}^{*}-\beta_{\text{in
}}\nabla_{\omega_{m-1}^{*}}\ell_{\mathrm{SSM}}\left(\mathcal{D}_{\text{target}};\omega_{m-1}^{*}\right),\end{split}$
(25)

with $\omega_{0}^{*}$ being the parameters initialized for the NSSMs.

It is evident from Fig. 3 that meta-learning expedites the training of NSSMs
by using data from similar systems. Moreover, Fig. 5 shows the tracking
performance of different algorithms. Here, we aim to track a circle (centered
at the origin with radius $2$) by running MPC on the NSSMs trained for $0$,
$100$, and $3000$ steps, respectively. From the first column of the figure, we
conclude that after the meta-training phase, both iMAML and MAML produce good
aggregated models that can capture dynamics of the target system to a certain
extent. The performance is further improved after the adaptation steps using
target data $\operatorname*{\mathcal{D}}_{\text{target}}$. However, since
supervised learning does not benefit from the similar systems, it performs
poorly compared with using meta-learning.

Another observation from comparing iMAML and MAML in Fig. 3 and Fig. 5 is that
iMAML tends to exhibit faster convergence and better tracking performance,
especially in the initial stages of meta-inference. The enhancement is
particularly striking when comparing iMAML and MAML tracking after 10 and 100
steps. This is because that iMAML exhibits superior performance with lower
loss during the meta-training phase, as evidenced by Fig. 2. Consequently, it
starts with better initial parameters, facilitating more effective adaptation
to the target system.

Figure 3: Comparison of the prediction performance on the target system.
Figure 4: The Pendulum Gym Environment.

(a) The iMAML-based MPC (Algorithm 2).

(b) The MAML-based MPC.

(c) The supervised learning-based MPC.

Figure 5: Comparison of the tracking performance on the target system, where
the three columns show the performance of different algorithms by using the
NSSMs adapted after $10$, $100$, and $3000$ steps, respectively. The lines and
shaded regions respectively represent the mean and standard deviation across
$10$ runs.

###  IV-B Pendulum

We also test the proposed algorithm on the benchmark Pendulum Gym Environment
[24]. The environment consists of a pendulum that is free to swing in a
vertical plane under the influence of gravity and is controlled by applying a
torque to the end point of the pendulum, see Fig. 4. The goal is to keep the
pendulum upright, with the least amount of applied torque. The state
$x=[\theta;\dot{\theta}]$ consists of the angle and angular velocity of the
pendulum, while the input is the torque $\tau$ and the output is $\theta$. We
collect the input and output data that are contaminated by random noises.

We construct the source systems by sampling pendulum masses of varying weights
and then evaluate on five target tasks involving pendulums of different mass.
Due to the space limits, we omit parameters in Algorithm 1 and 2. Fig. 6 shows
the prediction performance of NSSM along the meta-training phase. When this
phase converges, we apply the aggregated model on target systems. For all
target systems, the control error, defined as the deviation of the angle from
the center, i.e., $|\theta|$, decreases to below $10^{-4}$ after 100
adaptation steps. That means, our algorithms are robust to noise and
successfully stablize the pendulum by using few samples on target systems.
Therefore, meta learning is useful to pre-train a reasonable NSSM that can be
adapted to system changes and quickly finds proper controllers.

Figure 6: The prediction performance of NSSM on the pendulum.

##  V Conclusion

This paper considers the problem of reference tracking in an unknown nonlinear
system by using an NSSM. It introduces an iMAML-based MPC algorithm, which
comprises two phases. In the meta-training phase, data from source systems is
leveraged to pre-train the NSSM, while in the meta-inference phase, this model
is quickly adapted to the target system using only a small amount of data.
This approach addresses the issue of limited dataset availability for the
target system, enabling the development of a precise model that facilitates
optimal tracking even with a small dataset. Numerical examples demonstrate the
superior performance of the proposed algorithm compared to several baselines.
As a future work, we are planning experiments on physical benchmark problems
with a view towards transitioning to manufacturing processes.

## References

  * [1] F. Lin and R. D. Brandt, “An optimal control approach to robust control of robot manipulators,” _IEEE Transactions on Robotics and Automation_ , vol. 14, no. 1, pp. 69–77, 1998.
  * [2] H. Pan and M. Xin, “Nonlinear robust and optimal control of robot manipulators,” _Nonlinear Dynamics_ , vol. 76, pp. 237–254, 2014.
  * [3] D. L. Pepyne and C. G. Cassandras, “Optimal control of hybrid systems in manufacturing,” _Proceedings of the IEEE_ , vol. 88, no. 7, pp. 1108–1123, 2000.
  * [4] T. Dierks and S. Jagannathan, “Online optimal control of affine nonlinear discrete-time systems with unknown internal dynamics by using time-based policy update,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 23, no. 7, pp. 1118–1129, 2012.
  * [5] H. Modares, F. L. Lewis, W. Kang, and A. Davoudi, “Optimal synchronization of heterogeneous nonlinear systems with unknown dynamics,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 1, pp. 117–131, 2017.
  * [6] B. Zhao, D. Liu, and C. Luo, “Reinforcement learning-based optimal stabilization for unknown nonlinear systems subject to inputs with uncertain constraints,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 31, no. 10, pp. 4330–4340, 2019.
  * [7] A. Chakrabarty, G. Wichern, and C. R. Laughman, “Meta-learning of neural state-space models using data from similar systems,” _IFAC-PapersOnLine_ , vol. 56, no. 2, pp. 1490–1495, 2023.
  * [8] C. Legaard, T. Schranz, G. Schweiger, J. Drgoňa, B. Falay, C. Gomes, A. Iosifidis, M. Abkar, and P. Larsen, “Constructing neural network based models for simulating dynamical systems,” _ACM Computing Surveys_ , vol. 55, no. 11, pp. 1–34, 2023.
  * [9] D. Masti and A. Bemporad, “Learning nonlinear state–space models using autoencoders,” _Automatica_ , vol. 129, p. 109666, 2021.
  * [10] J. A. Suykens, B. L. De Moor, and J. Vandewalle, “Nonlinear system identification using neural state space models, applicable to robust control design,” _International Journal of Control_ , vol. 62, no. 1, pp. 129–152, 1995.
  * [11] B. O. Koopman and J. v. Neumann, “Dynamical systems of continuous spectra,” _Proceedings of the National Academy of Sciences_ , vol. 18, no. 3, pp. 255–263, 1932.
  * [12] M. Korda and I. Mezić, “Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control,” _Automatica_ , vol. 93, pp. 149–160, 2018.
  * [13] L. Xin, L. Ye, G. Chiu, and S. Sundaram, “Identifying the dynamics of a system by leveraging data from similar systems,” in _2022 American Control Conference (ACC)_. IEEE, 2022, pp. 818–824.
  * [14] A. Chakrabarty, “Optimizing closed-loop performance with data from similar systems: A bayesian meta-learning approach,” in _2022 IEEE 61st Conference on Decision and Control (CDC)_. IEEE, 2022, pp. 130–136.
  * [15] L. Li, C. De Persis, P. Tesi, and N. Monshizadeh, “Data-based transfer stabilization in linear systems,” _IEEE Transactions on Automatic Control_ , 2023.
  * [16] V. Dumoulin, N. Houlsby, U. Evci, X. Zhai, R. Goroshin, S. Gelly, and H. Larochelle, “Comparing transfer and meta learning approaches on a unified few-shot classification benchmark,” _arXiv preprint arXiv:2104.02638_ , 2021.
  * [17] L. F. Toso, D. Zhan, J. Anderson, and H. Wang, “Meta-learning linear quadratic regulators: A policy gradient maml approach for the model-free lqr,” _arXiv preprint arXiv:2401.14534_ , 2024.
  * [18] S. Zhan, G. Wichern, C. Laughman, A. Chong, and A. Chakrabarty, “Calibrating building simulation models using multi-source datasets and meta-learned bayesian optimization,” _Energy and Buildings_ , vol. 270, p. 112278, 2022.
  * [19] R. Kaushik, T. Anne, and J.-B. Mouret, “Fast online adaptation in robotics through meta-learning embeddings of simulated priors,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 5269–5276.
  * [20] X. Song, Y. Yang, K. Choromanski, K. Caluwaerts, W. Gao, C. Finn, and J. Tan, “Rapidly adaptable legged robots via evolutionary meta-learning,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 3769–3776.
  * [21] A. Rajeswaran, C. Finn, S. M. Kakade, and S. Levine, “Meta-learning with implicit gradients,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [22] C. E. Garcia, D. M. Prett, and M. Morari, “Model predictive control: Theory and practice—a survey,” _Automatica_ , vol. 25, no. 3, pp. 335–348, 1989.
  * [23] C. Finn, P. Abbeel, and S. Levine, “Model-agnostic meta-learning for fast adaptation of deep networks,” in _International conference on machine learning_. PMLR, 2017, pp. 1126–1135.
  * [24] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, “Openai gym,” _arXiv preprint arXiv:1606.01540_ , 2016.