# Absolute Properties of An Overcontact Binary HH Boo

H.A. Dal<EMAIL_ADDRESS>E. Sipahi  Ege University, Science Faculty,
Department of Astronomy and Space Sciences, 35100 Bornova, İzmir, Turkey

###### Abstract

We obtained multi-colour light curves of HH Boo. We analysed the orbital
period variation of the system. The analysis indicated that there is possible
mass transfer from the second component to the primary or mass loss with
$-5.04\times 10^{-7}$ $M_{\odot}$ per year. Re-analysing the available radial
velocity curve, we analysed the light curves. The inclination ($i$) of the
system was found to be 69∘.71$\pm$0∘.16, while the semi-major axis ($a$) was
computed as 2.246$\pm$0.064 $R_{\odot}$. The mass of the primary component was
found to be 0.92$\pm$0.08 $M_{\odot}$, while it was obtained as 0.58$\pm$0.06
$M_{\odot}$ for the secondary component. The radius of the primary component
was computed as 0.98$\pm$0.03 $R_{\odot}$, while it was computed as
0.80$\pm$0.02 $R_{\odot}$ for the secondary component. We demonstrated that HH
Boo should be a member of the A-type subclass of W UMa binaries.

###### keywords:

techniques: photometric — (stars:) binaries: eclipsing — stars: late-type —
stars: individual: (HH Boo)

††journal: New Astronomy

##  1 Introduction

HH Boo (GSC 03472-00641) is classified as an eclipsing binary of W UMa type (a
contact binary) in the SIMBAD Database. For the first time, the system was
listed in the TYCHO-2 Catalogue by Høg et al. (2000), who gave $V=11^{m}.32$
and $(B-V)=0^{m}.45$ for the system. Its variability nature was found by
Maciejewski et al. (2003), who obtained the first light curve and gave the
light elements as following: $T_{0}=2452764.50965$ and $P=0^{d}.318618$.
Considering the spectra of the system, Maciejewski et al. (2003) indicated
that the system should be from the spectral type G5III. Then, Maciejewski &
Ligeza (2004) obtained the radial velocity curve of the system. They found
that mass ratio is 0.633$\pm$0.042, while $M_{1}sin^{3}i=0.78\pm 0.08$
$M_{\odot}$ and $M_{2}sin^{3}i=0.49\pm 0.05$ $M_{\odot}$. Although Maciejewski
et al. (2003) have obtained a light curve, there is no light curve analysis
given in the literature. However, many minima times have been obtained along
the years. Most of them are listed in the modern database of O-C Gateway
(Paschke & Brát, 2006), while some others are given by Maciejewski & Karska
(2004), Maciejewski & Niedzielski (2004), Brát et al. (2008), Hubscher et al.
(2009), Hubscher et al. (2010), Walter (2010).

According to Ammons et al. (2006), who determined new temperatures and
metallicities for more than 100,000 FGK dwarfs, the temperature of the system
is 5699 K, while $[Fe/H]$ was found to be -0.57. The distance is 54 pc.

In the literature, there are several systems similar to HH Boo (Essam et al.,
2010). In this study, we obtained multi-colour light curves. Analysing the
orbital period variation, we adjusted the light elements of the system and
derived the mass transfer rate between the components. Then, re-analysing the
radial velocity obtained by Maciejewski & Ligeza (2004), we analysed the
multi-colour light curve using the parameters such as mass ratio and semi-
major axis.

##  2 Observations

Observations of the system were acquired with a thermoelectrically cooled ALTA
$U+47$ 1024$\times$1024 pixel CCD camera attached to a 35 cm - Schmidt -
Cassegrains - type MEADE telescope at Ege University Observatory. The
observations were continued in BVR bands in the observing season 2011. Some
basic parameters of program stars are listed in Table 1, in which the
brightness and colours were taken from the General Catalogue of Variable Stars
(Kholopov et al., 1988).

Although the program and comparison stars are very close on the sky,
differential atmospheric extinction corrections were applied. The atmospheric
extinction coefficients were obtained from observations of the comparison
stars on each night. Heliocentric corrections were also applied to the times
of the observations. The mean averages of the standard deviations are 0m.023,
0m.011, and 0m.010 for observations acquired in the BVR bands, respectively.
To compute the standard deviations of observations, we used the standard
deviations of the reduced differential magnitudes in the sense comparisons
(GSC 03472-00043) minus check (GSC 03472-01201) stars for each night. There
was no variation observed in the standard brightness comparison stars.

##  3 Analysis of Orbital Period Variation

There are several minima times of the system in the literature (Paschke &
Brát, 2006; Hubscher et al., 2010). In the analysis, we used nor visual
observations neither the minima times with large error. We also obtained new
15 minima times. All the minima times used in this study are listed in Table
2. The standard deviation of each one obtained in this study is given in the
brackets near its related digits in Table 2.

Using the regression calculation, investigations demonstrated that the main
variation could be described by a downward parabolic curve, which must be
caused by possible mass transfer from the second component to the primary or
mass loss. Therefore, the main period variation was represented by the
quadratic light elements, which are given by Equation (1).

$Min~{}I~{}(Hel.)~{}=~{}24~{}54912.3236(2)~{}+~{}0^{d}.318666(1)~{}\times~{}E-1.3349(1)\times
10^{-11}~{}\times~{}E^{2}$ (1)

where the standard deviations of each coefficient and each constant are given
in the brackets near their related digits.

The period variation and its quadratic fit are shown in Figure 1. In the
analyses, the weighted sum of the squared residuals, $\Sigma w(O-C)^{2}$, was
found to be 0.000038 $day^{2}$. Considering the quadratic term ($Q$), the
parameter of the period variation ($dP/dt$) was found to be -9.60$\times
10^{-7}$ $yr^{-1}$. The relation between the parameter of the period variation
($dP/dt$) and the rate of mass transfer was determined with Equation (2) by
Pringle & Wade (1985):

$\dot{m}~{}=~{}\frac{1}{3}~{}\times~{}\frac{M_{1}~{}\times~{}M_{2}}{M_{1}~{}-~{}M_{2}}~{}\times~{}\frac{\dot{P}}{P}$
(2)

where $\dot{m}$ is mass transfer rate per year. Using this equation, the mass
transferring from the secondary component to the primary was found to be
$-5.04\times 10^{-7}$ $M_{\odot}~{}yr^{-1}$.

##  4  Radial Velocity and Light Curve Analyses

We took the radial velocity curve of the system from Maciejewski & Ligeza
(2004). In order to adjust their radial velocity solution, we used the
Spectroscopic Binary Solver software (Johnson, 2004) and re-analysed the
radial velocity curve. In the analysis, we used the orbital period adjusted by
the $O-C$ analyses in this study, and we fixed it. The results are listed in
Table 3. The solution generally gave the same values obtained by Maciejewski &
Ligeza (2004) with the little differences. The theoretical radial velocity
curves derived by the Spectroscopic Binary Solver software are shown in Figure
2.

HH Boo is classified as an eclipsing binary of W UMa type in the SIMBAD
Database, and it is more likely a contact binary system. In this respect, we
analysed the light curves obtained in the BVR bands with using the PHOEBE
V.0.31a software (Prša & Zwitter, 2005), which is used in the version 2003 of
the Wilson-Devinney Code (Wilson & Devinney, 1971; Wilson, 1990). We tried to
analyse the light curves with different modes, such as the ”overcontact binary
not in thermal contact”, ”semi-detached system with the primary component
filling its Roche-Lobe”, and ”double contact binary” modes. The initial
analyses demonstrated that an astrophysical acceptable result can be obtained
if the analysis is carried out in the ”overcontact binary not in thermal
contact” mode, while no acceptable results in the astrophysical sense could be
statistically obtained in all the others modes.

According to the radial velocity analysis, the mass ratio of the system should
be 0.632$\pm$0.038. Ammons et al. (2006) determined the temperature of the
system as 5699 K. Thus, the temperature of the primary component was fixed to
5699 K in the analyses, and the temperature of the secondary was taken as a
free parameter. Considering the spectral type corresponding to this
temperature, the albedos ($A_{1}$ and $A_{2}$) and the gravity darkening
coefficients ($g_{1}$ and $g_{2}$) of the components were adopted for the
stars with the convective envelopes (Lucy, 1967; Rucinski, 1969). The non-
linear limb-darkening coefficients ($x_{1}$ and $x_{2}$) of the components
were taken from Van Hamme (1993). In the analyses, the fractional luminosity
($L_{1}$) of the primary component and the inclination ($i$) of the system
were also taken as the adjustable free parameters.

The parameters derived from the analyses are listed in Table 4, while the
synthetic light curves are shown in Figure 3. In addition, using the
parameters obtained from the light curve analysis, we also derived the Roche
geometry of the system that is shown in Figure 4.

Using the radial velocity curves, the mass ratio of the system ($q$) was
obtained as 0.632$\pm$0.038, and the semi-major axis ($a$) was found to be
2.246$\pm$0.064 $R_{\odot}$. Considering both the radial velocity curve
solution and the inclination ($i$) of the system found from the light curve
analysis, the masses were found to be 0.92$\pm$0.08 $M_{\odot}$ for the
primary component and 0.58$\pm$0.06 $M_{\odot}$ for the secondary component.
Considering the semi-major axis, the radius of the primary component was
computed as 0.98$\pm$0.03 $R_{\odot}$, while it was computed as 0.80$\pm$0.02
$R_{\odot}$ for the secondary component. In addition, the luminosity of the
primary component was computed as 0.91$\pm$0.08 $L_{\odot}$, and it was
computed as 0.47$\pm$0.03 $L_{\odot}$ for the secondary component.

In order to test whether the absolute parameters are generally acceptable in
the astrophysical sense, or not, we compared the components with other systems
in the mass-radius ($M-R$), mass-luminosity ($M-L$), and luminosity-effective
temperature ($T_{eff}-L$) planes. All the comparisons are shown in Figure 5.
In the figure, the lines represent the ZAMS theoretical model developed for
the stars with $Z=0.02$ by Girardi et al. (2000), while dashed lines represent
the TAMS theoretical model. The filled circles represent the primary
components, while the open circles represents the secondary ones. The
components of HH Boo are located together with some samples of its analogues,
such as YY CrB, DN Boo, CK Boo, $\epsilon$ CrA, FG Hya, TV Mus, AW UMa, GR
Vir, V776 Cas. The sample systems were taken from Essam et al. (2010), and
they are shown in purple colour, while HH Boo is shown in black colour in
Figure 5.

##  5 Discussion

In this study, we tried to determine the nature of an eclipsing binary system
HH Boo. The analysis of orbital period variation indicates possible mass
transfer from the second component to the primary and/or possible mass loss
from the system. The size of transferring or loosing mass was computed as
$-5.04\times 10^{-7}$ $M_{\odot}$ per year, which causes a decreasing in the
orbital period due to the angular momentum loss. In fact, the parameter of the
period variation was found to be -9.60$\times 10^{-7}$ $yr^{-1}$. We adjusted
the orbital period as $0^{d}.318666$. It is well known that several binaries
of W UMa type such as YY CrB (Essam et al., 2010), BS Cas (Yang et al., 2008),
VZ Tri (Yang, 2010), exhibit large the period variation due to the large mass
transfer.

We also re-analysed available radial velocity curve. Our solution seems to be
a little bit different according to the results found by Maciejewski & Ligeza
(2004). The differences should be caused due to the orbital period used in
this study. Here, the adjusted period is a bit different from the one used by
Maciejewski & Ligeza (2004).

For the first time in the literature, we analysed light curves of the system.
The inclination ($i$) of the system was found to be 69∘.71$\pm$0∘.16, while
the temperature of the secondary component was found to be 5352$\pm$15 K from
the analysis. The fractional radii were found to be $r_{1}=0.435\pm 0.001$ for
the primary component and $r_{2}=0.354\pm 0.001$ for the secondary one. In
this case, the sum of fractional radii was computed as $r_{1}+r_{2}\simeq
0.80$. Thus, HH Boo seems to be in agreement with Kopal (1956)’s criteria for
overcontact systems. The $O-C$ analysis indicates that the orbital period as
$0^{d}.318666$. In addition, the temperature of the primary component is 5699
K, while the secondary one is 5352 K. Although some W UMa type binaries have
components with some different surface temperature, they generally have the
same surface temperature. Here, the primary component of HH Boo is a little
bit hotter than the secondary one. The $O-C$ analysis indicates possible mass
transfer from the second component to the primary. This case should be the
reason of the hotter primary component. Considering some characteristics of
the system such as the short orbital period, small mass ratio, hotter primary
component and ect., HH Boo seems to be in agreement with the members of the
A-type subclass of W UMa binaries (Berdyugina, 2005; Rucinski, 1985).

As it is seen from Figure 5, comparing HH Boo with its analogues in some
planes, such as $M-R$, $M-L$, and $T_{eff}-L$ planes, demonstrated that the
components of the system are in agreement with their analogues. On the other
hand, both components are seen closer to each other in the figures. If the
result obtained from the $O-C$ analysis is taken into account, it is possible
that the secondary component will continue to lose its mass, while the primary
will collect. Therefore, both components will have been separated from each
other in the planes shown in Figure 5. However, they will get closer to each
other in their orbits as it is in the models of Rocha-Pinto et al. (2002).

## Acknowledgments

The author acknowledges the generous observing time awarded to the Ege
University Observatory. We also thank the referee for useful comments that
have contributed to the improvement of the paper.

## References

  * Ammons et al. (2006) Ammons, S.M., Robinson, S.E., Strader, J., Laughlin, G., Fischer, D., Wolf, A., 2006, ApJ, 638, 1004
  * Berdyugina (2005) Berdyugina, S.V., 2005, LRSP, 2, 8
  * Brát et al. (2008) Brát, L., Šmelcer, L., Kuèáková, H., Ehrenberger, R., Kocián, R., Lomoz, F., Urbanèok, L., Svoboda, P., Trnka, J., Marek, P., and 5 coauthors, 2008, OEJV, 94, 1
  * Essam et al. (2010) Essam, A., Saad, S.M., Nouh, M.I., Dumitrescu, A., El-Khateeb, M.M., Haroon, A., 2010, NewA, 15, 227
  * Girardi et al. (2000) Girardi, L., Bressan, A., Bertelli, G., Chiosi, C., 2000, A&AS 141, 371
  * Høg et al. (2000) Høg, E., Fabricius, C., Makarov, V.V., Urban, S., Corbin, T., Wycoff, G., Bastian, U., Schwekendiek, P., Wicenec, A., 2000, A&A, 355L, 27
  * Hubscher et al. (2010) Hubscher, J., Lehmann, P.B., Monninger, G., Steinbach, H., Walter, F., 2010, IBVS, 5918, 1
  * Hubscher et al. (2009) Hubscher, J., Steinbach, H.M., Walter, F., 2009, IBVS, 5874, 1
  * Johnson (2004) Johnson, D.O., 2004, JAD, 10, 3
  * Kholopov et al. (1988) Kholopov, P.N., Samus, N.N., Frolov, M.S., Goranskij, V.P., Gorynya, N.A., Karitskaya, E.A., Kazarovets, E.V., Kireeva, N.N., et al., 1985-1988, General Catalogue of Variable Stars, Vol. I III, 4th ed., Nauka, Moscow
  * Kopal (1956) Kopal, Z., 1956. AnAp, 19, 298
  * Lucy (1967) Lucy, L.B., 1967, Z. Astrophys, 65, 89
  * Maciejewski et al. (2003) Maciejewski, G., Czart, K., Niedzielski, A., Karska, A., 2003, IBVS, 5431, 1
  * Maciejewski & Karska (2004) Maciejewski, G., Karska, A., 2004, IBVS, 5494, 1
  * Maciejewski & Ligeza (2004) Maciejewski, G., Ligeza, P., 2004, IBVS, 5504, 1
  * Maciejewski & Niedzielski (2004) Maciejewski, G., Niedzielski, A., 2004, BaltA, 13, 700
  * Paschke & Brát (2006) Paschke, A., Brát, L., 2006, OEJV, 23, 13
  * Pringle & Wade (1985) Pringle, J.E. and Wade, R.A., 1985, ”Interacting Binary Stars”, Cambridge University Press., New York, p.17
  * Prša & Zwitter (2005) Prša, A., Zwitter, T., 2005, ApJ, 628, 426
  * Rocha-Pinto et al. (2002) Rocha-Pinto, H.J., Castilho, B.V., Maciel,W.J., 2002, A&A, 384, 912
  * Rucinski (1969) Rucinski, S.M., 1969, AcA, 19, 245
  * Rucinski (1985) Rucinski, S., 1985. In: Pringle, J.E., Wade, R.A. (Eds.), Interacting binary stars. Cambridge University Press, Cambridge. p. 1
  * Yang et al. (2008) Yang, Y.-G., Wei, J.-Y., He, J.-J., 2008, AJ, 136, 594
  * Yang (2010) Yang, Y.-G., 2010, Ap&SS, 326, 125
  * Van Hamme (1993) Van Hamme, W., 1993, AJ, 106, 2096
  * Walter (2010) Walter, F., 2010, BAVSR, 59, 139
  * Wilson (1990) Wilson, R.E., 1990, ApJ, 356, 613
  * Wilson & Devinney (1971) Wilson, R.E., Devinney, E.J., 1971, ApJ, 166, 605

Figure 1: The O-C variation of HH Boo. Figure 2: Radial velocity curve of HH
Boo. Open circles represent the observations of the primary, while filled
circles represent the secondary components. Solid curves are the theoretical
radial velocity curves derived by the Spectroscopic Binary Solver software.
Figure 3: The BVR light curves of HH Boo and the synthetic solutions for the
observations in each band. Figure 4: The Roche geometry of HH Boo. Figure 5:
The places of the components of HH Boo in the planes of (upper panel) the
mass-radius ($M-R$), (middle panel) mass-luminosity ($M-L$), and (bottom
panel) luminosity-effective temperature ($T_{eff}-L$). In the panels, the
continuous and dashed lines represent the ZAMS and TAMS theoretical models
developed by Girardi et al. (2000), respectively. The filled circles represent
the primary components, while the open circles represent the secondary ones.
The dark circles represent the HH Boo components, while the purple coloured
circles represent the components of other contact binaries.

Table 1: Basic parameters for the observed stars. Star | Alpha (J2000) | Delta (J2000) | V | B-V
---|---|---|---|---
| (h m s) | (∘ ′ ′′) | (mag) | (mag)
HH Boo | 14 21 44.06 | 46 41 59.40 | 11.272 | 0.654
GSC 03472-00043 | 14 21 46.59 | 46 43 49.40 | 12.906 | 0.595
GSC 03472-01201 | 14 21 37.30 | 46 46 18.90 | 13.701 | 0.662
Table 2: The minima times and ($O-C$) residuals (In the first column, the standard deviations of obtained minima times are given in the brackets near themselves). O | E | $(O-C)_{II}$ | Type | Method | REF
---|---|---|---|---|---
52749.5328 | -9420.0 | -0.00179 | I | V | Paschke & Brát (2006)
52764.5096 | -9373.0 | -0.00230 | I | V | Paschke & Brát (2006)
54148.6406 | -5029.5 | 0.00192 | II | Ir | Paschke & Brát (2006)
54513.5136 | -3884.5 | 0.00209 | II | Ir | Hubscher et al. (2010)
54599.3935 | -3615.0 | 0.00146 | I | CCD+R | Paschke & Brát (2006)
54912.3233 | -2633.0 | 0.00100 | I | Ir | Paschke & Brát (2006)
54912.4835 | -2632.5 | 0.00187 | II | Ir | Paschke & Brát (2006)
54937.3392 | -2554.5 | 0.00160 | II | Ir | Paschke & Brát (2006)
54958.6887 | -2487.5 | 0.00046 | II | V | Paschke & Brát (2006)
54958.8490 | -2487.0 | 0.00143 | I | V | Paschke & Brát (2006)
55644.9336 | -334.0 | -0.00237 | I | V | Paschke & Brát (2006)
55680.7853 | -221.5 | -0.00062 | II | V | Paschke & Brát (2006)
55716.4767(3) | -109.5 | 0.00013 | II | R | This Study
55716.4768(5) | -109.5 | 0.00023 | II | B | This Study
55716.4768(2) | -109.5 | 0.00025 | II | V | This Study
55743.4026(3) | -25.0 | -0.00128 | I | R | This Study
55743.4032(4) | -25.0 | -0.00060 | I | V | This Study
55743.4038(2) | -25.0 | -0.00008 | I | B | This Study
55747.3863(5) | -12.5 | -0.00083 | II | R | This Study
55747.3871(3) | -12.5 | -0.00003 | II | V | This Study
55747.3874(4) | -12.5 | 0.00022 | II | B | This Study
55751.3696(3) | 0.0 | -0.00087 | I | V | This Study
55751.3700(3) | 0.0 | -0.00051 | I | R | This Study
55751.3700(6) | 0.0 | -0.00050 | I | B | This Study
55751.3701(2) | 0.0 | -0.00037 | I | R | This Study
55751.3702(3) | 0.0 | -0.00027 | I | V | This Study
55751.3702(4) | 0.0 | -0.00026 | I | B | This Study
Table 3: The results of the analysis of Radial Velocity curve. Parameter | Value | Error
---|---|---
Long. of Periastron1 (∘) | 238.669 | ${}^{+49.884}_{-72.41}$
Long. of Periastron2 (∘) | 58.669 | ${}^{+49.884}_{-72.41}$
Eccentricity ($e$) | 0.048 | ${}^{+0.043}_{-0.042}$
Semi-Amplitude1 ($kms^{-1}$) | 129.647 | ${}^{+6.4243}_{-6.2471}$
Semi-Amplitude2 ($kms^{-1}$) | 205.078 | ${}^{+7.0675}_{-6.7157}$
Systemic Velocity ($kms^{-1}$) | 2.9209 | ${}^{+4.0328}_{-4.0429}$
$a_{1}sini$ ($km$) | $5.67\times 10^{5}$ | $\pm 3.68\times 10^{4}$
$a_{2}sini$ ($km$) | $8.98\times 10^{5}$ | $\pm 4.86\times 10^{4}$
$m_{1}sin^{3}i$ ($M_{\odot}$) | $7.56\times 10^{-1}$ | $\pm 9.99\times 10^{-2}$
$m_{2}sin^{3}i$ ($M_{\odot}$) | $4.78\times 10^{-1}$ | $\pm 6.54\times 10^{-2}$
Table 4: The parameters of components obtained from the light curve analysis. Parameter | Value
---|---
$i$ (∘) | 69.71$\pm$0.16
$T_{1}$ (K) | 5699
$T_{2}$ (K) | 5352$\pm$15
$\Omega_{1}$ | 3.0492
$\Omega_{2}$ | 3.0402$\pm$0.0051
$L_{1}/L_{T}$ (B) | 0.691$\pm$0.046
$L_{1}/L_{T}$ (V) | 0.669$\pm$0.037
$L_{1}/L_{T}$ (R) | 0.657$\pm$0.031
$g_{1}$, $g_{2}$ | 0.32, 0.32
$A_{1}$, $A_{2}$ | 0.50, 0.50
$x_{1,bol}$, $x_{2,bol}$ | 0.619, 0.619
$x_{1,B}$, $x_{2,B}$ | 0.765, 0.765
$x_{1,V}$, $x_{2,V}$ | 0.732, 0.732
$x_{1,R}$, $x_{2,R}$ | 0.655, 0.655
$<r_{1}>$ | 0.435$\pm$0.001
$<r_{2}>$ | 0.354$\pm$0.001


# Causal Repair of Learning-Enabled
Cyber-Physical Systems

Pengyuan Lu  Computer and Information Science
University of Pennsylvania
Philadelphia, PA, USA
<EMAIL_ADDRESS>Ivan Ruchkin  Electrical and Computer Engineering
University of Florida
Gainesville, FL, USA
<EMAIL_ADDRESS>Matthew Cleaveland  Oleg Sokolsky  Computer and
Information Science
University of Pennsylvania
Philadelphia, PA, USA
<EMAIL_ADDRESS>Computer and Information Science
University of Pennsylvania
Philadelphia, PA, USA
<EMAIL_ADDRESS>Insup Lee  Computer and Information Science
University of Pennsylvania
Philadelphia, PA, USA
<EMAIL_ADDRESS>
###### Abstract

Models of actual causality leverage domain knowledge to generate convincing
diagnoses of events that caused an outcome. It is promising to apply these
models to diagnose and repair run-time property violations in cyber-physical
systems (CPS) with learning-enabled components (LEC). However, given the high
diversity and complexity of LECs, it is challenging to encode domain knowledge
(e.g., the CPS dynamics) in a scalable actual causality model that could
generate useful repair suggestions. In this paper, we focus causal diagnosis
on the input/output behaviors of LECs. Specifically, we aim to identify which
subset of I/O behaviors of the LEC is an actual cause for a property
violation. An important by-product is a counterfactual version of the LEC that
repairs the run-time property by fixing the identified problematic behaviors.
Based on this insights, we design a two-step diagnostic pipeline: (1)
construct and Halpern-Pearl causality model that reflects the dependency of
property outcome on the component’s I/O behaviors, and (2) perform a search
for an actual cause and corresponding repair on the model. We prove that our
pipeline has the following guarantee: if an actual cause is found, the system
is guaranteed to be repaired; otherwise, we have high probabilistic confidence
that the LEC under analysis did not cause the property violation. We
demonstrate that our approach successfully repairs learned controllers on a
standard OpenAI Gym benchmark.

###### Index Terms:

actual causality, control policy repair, cyber-physical system.

##  I Introduction

When a person’s leg hurts, they seek detailed diagnosis from a physician
regarding the pain’s cause, which would lead to an effective intervention as a
“repair”. Similarly, when a closed-loop cyber-physical system (CPS) violates a
desirable property at run time, the violation needs diagnosis — a procedure
that identifies the cause for the violation, and by fixing the identified
cause, the CPS can be repaired. Traditionally, researchers conduct this kind
of analysis from statistical inference on observations [1, 2] . However,
statistical diagnosis is prone to mistaking correlation for causation: when a
student always wears a green jacket and fails several exams, such algorithms
are likely to conclude it is the green jacket’s fault due to the perfect
correlation. Therefore, in this paper, we focus on stronger causal reasoning
and repair on CPS failures.

Researchers have explored the concept of actual causality that leverages
domain knowledge to produce well-defined and convincing causal explanations.
Informally, an actual cause for an outcome is a minimal set of variable
assignments, which represent an event, that changes the outcome if assigned
some counterfactual values. Finding an actual cause requires the construction
of an actual causality model, such as a Halpern-Pearl model, which rigorously
defines actual causes for events and precisely assigns responsibility and
blame [3, 4, 5]. With these models, we can encode common knowledge that, for
instance, the student’s bad grade can be either due to not studying hard
enough or misunderstanding some concepts in class.

Although actual causality is promising for analyzing and fixing CPS failures,
causal analysis and repair have been complicated by the growing popularity of
learning-enabled components (LECs) [6, 7, 8]. First, LECs usually consist of
numerous internal continuous parameters, such as weights and biases in deep
neural networks. Second, LECs take a large diversity of forms, from basic
statistical models such as linear regressors and support vector machines to
complex deep architectures, lacking a shared struture of a parameter template.
Furthermore, sometimes the internal structures of LECs are black-boxes
protected as intellectual properties, and are invisible to testing engineers.
Under these scenarios, it is hard to build a causal model based on the
internal information flow of these components. Related to diagnosis and repair
are the efforts in explainable AI [9, 10] and formal methods [11, 12], where
researchers build frameworks to explain behaviors of learned agents. However,
to our knowledge, these strands of work have yet to connect to actual
causality.

Since the internal parameters of LECs are hard to analyze, we take a step back
and analyze the LEC I/O behavior instead. Our intention in this paper is to
leverage actual causality to efficiently identify the granular I/O behaviors
of a suspected LEC for a run-time property violation. We intend to identify
these behaviors in the form of fine-grained input-output mappings, so that a
precise repair can be made for CPS to satisfy the desired property.

To achieve this goal, we first construct an Halpern-Pearl (HP) causality model
to encode the dependency of property outcome on the suspected LEC’s behaviors.
Then, we design an algorithm to search for an actual cause and corresponding
repair on this HP model. This algorithm either concludes that an actual cause
does not exist within the LEC’s behaviors with high confidence , or outputs
the set of behaviors that are causing the violation, along with a
counterfactual repair. Experiments on the mountain car, a standard OpenAI Gym
benchmark [13], demonstrates the capability of our method.

To summarize, our contributions are:

  1. 1.

A novel use of HP causality model to identify problematic LEC I/O behaviors
that cause a run-time property violation, and produce repair suggestions,
without touching the component’s internal information flow, and

  2. 2.

Experiments on an OpenAI Gym benchmark that show our solution’s utility for
learning-enabled CPS repair.

##  II Background and Related Work

###  II-A Learning-enabled Components (LEC)

Learning-enabled components (LECs) are functional components of larger systems
that are learned from data, either offline or online. One example is the
perception unit, with neural networks able to complete difficult vision tasks
that generally cannot be accomplished by traditional first-principles
algorithms, such as tracking moving objects at 100 fps [6]. Moreover, deep
reinforcement learning has provided useful control policies on various tasks
as a substitute for conventional controllers [14, 15]. Due to LECs’
performance, the state-of-the-art in CPS has a growing enthusiasm of these
statistically generated agents, with automated design tools built for CPS with
LECs [7], as well as analysis on their assurance on run-time properties such
as safety [8]. Consequently, we need to consider the presence of LECs when
diagnosing faults in CPS.

###  II-B Repair

A repair is a procedure to change or replace parts of a system to achieve a
desirable performance, of which the system previously fell short. This term
has been widely used in traditional embedded systems. For example, researchers
have studied repair on hardware components such as DRAM [16, 17] or noisy
sensors [18]. Also, repair has been performed on software; for example,
assembly program can be transformed to decrease resource consumption below a
threshold [19].

In modern CPS with LEC, the topic of repair interests many researchers.
Learned agents, such as deep neural networks, may lead the system to unsafe
states [20] or unplanned paths [21]. Therefore, repair is necessary to recover
desirable behaviors. Repair of neural networks is an active area of research
[22, 23]. For instance, network parameters can be repaired by search
algorithms [24] or constraint solvers [25]. More recently, researchers have
studied provable repairs on deep neural networks, with guaranteed satisfaction
of a given property after the repair [26, 27, 28].

Compared to the existing work on repair, we focus on the causal relationships
between the repairable elements and the execution outcome. That is, the neural
network parameters to be repaired should be the ones that caused the system’s
failure.

###  II-C Actual Causality and Halpern-Pearl Models

Below we rephrase Halpern and Pearl’s definition of actual causality [3, 4].
An actual causality model, or a Halpern-Pearl (HP) model is a recursive
structure, i.e., a directed acyclic graph, such that every node represents
either an exogenous variable, whose value is determined by factors outside
this model, or an endogenous variable, whose value is determined by other
variables in this model. The edges represent dependencies: an exogenous node
has only outgoing edges but no incoming edges, while an endogenous node may
have both. Every endogenous node is equipped with a function, which defines
how the node’s value is computed from other nodes. In other words, this
function defines the incoming edges to the node. Formally,

  1. 1.

An _HP model_ is a tuple
$\mathcal{M}=(\mathcal{U}_{\text{endo}},\mathcal{U}_{\text{exo}},\mathcal{V},\mathcal{E})$,
where $\mathcal{U}_{\text{endo}}$ and $\mathcal{U}_{\text{exo}}$ are finite
sets of endogenous and exogenous variables, respectively, and for each
variable $u\in\mathcal{U}_{\text{endo}}\cup\mathcal{U}_{\text{exo}}$,
$\mathcal{V}(u)$ defines a non-empty and potentially infinite set of values
that $u$ can take.

  2. 2.

$\mathcal{E}$ is a set of edges, associated with dependency equations, that
defines how the value of each endogenous node $u$ is computed. I.e., for each
$e_{u}\in\mathcal{E}$,

$e_{u}:\prod_{u^{\prime}\in\mathcal{U}_{\text{endo}}\setminus\\{u\\}}\mathcal{V}(u^{\prime})\times\prod_{u^{\prime\prime}\in\mathcal{U}_{\text{exo}}}\mathcal{V}(u^{\prime\prime})\mapsto\mathcal{V}(u).$

To define an actual cause in an HP model $\mathcal{M}$, we introduce the
following notation:

  1. 1.

An assignment of a variable
$u\in\mathcal{U}_{\text{endo}}\cup\mathcal{U}_{\text{exo}}$ is denoted as
$u:=v$, for some $v\in\mathcal{V}(u)$. Assignments of multiple variables are
denoted in vector form $\mathbf{u}:=\mathbf{v}$, with
$\mathbf{u}=[u_{1},u_{2},\dots]$, $\mathbf{v}=[v_{1},v_{2},\dots]$. This
assignment means a conjunction $(u_{1}:=v_{1})\land(u_{2}:=v_{2})\land\dots$.

  2. 2.

A property $\psi$ is a Boolean function of endogenous variable assignments,
e.g., $\psi=(u_{1}:=v_{1})\land(u_{2}:=v_{2})\lor(u_{3}:=v_{3})$. The
satisfaction relation $(\mathcal{M},\mathbf{v}_{\text{cxt}})\models\psi$
denotes that a property holds on $\mathcal{M}$ given exogenous nodes assigned
with $\mathbf{v}_{\text{cxt}}$.

  3. 3.

A counterfactual $\mathbf{v}^{\prime}$ (with respect to the “factual”
$\mathbf{v}$) is an alternative value assignment on some endogenous variables.
The values of dependent nodes in $\mathbf{v}^{\prime}$ may be different from
$\mathbf{v}$ in accordance with to $\mathcal{E}$. A property $\psi$ holds on a
counterfactual that replaces the factual values $\mathbf{v}$ with
$\mathbf{v}^{\prime}$ on variables $\mathbf{u}$ is denoted as
$(\mathcal{M},\mathbf{v}_{\text{cxt}})\models[\mathbf{u}\leftarrow\mathbf{v}^{\prime}]\psi$,
or simply $[\mathbf{u}\leftarrow\mathbf{v}^{\prime}]\psi$.

Then, in the above terms, the Halpern and Pearl definition of an actual cause
is phrased as follows.

######  Definition 1 (Actual Cause).

On an HP model $\mathcal{M}$ with exogenous node values
$\mathbf{v}_{\text{cxt}}$, the assignment on a set of endogenous variables
$\mathbf{u}:=\mathbf{v}$ is an actual cause of $\psi$ iff the following
conditions hold:

  1. 1.

AC1:
$(\mathcal{M},\mathbf{v}_{\text{cxt}})\models(\mathbf{u}:=\mathbf{v})\land\psi$

  2. 2.

AC2: $\exists$ partition
$\mathcal{U}_{\text{endo}}=\mathbf{u}\cup\mathbf{u}_{1}\cup\mathbf{u}_{2}$.
Denote the factual values of $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ as
$\mathbf{v}_{1}$ and $\mathbf{v}_{2}$, respectively. Then,
$\exists\mathbf{v}^{\prime},\mathbf{v}_{1}^{\prime}$, such that
(a)
$[\mathbf{u}\leftarrow\mathbf{v}^{\prime},\mathbf{u}_{1}\leftarrow\mathbf{v}_{1}^{\prime}]\neg\psi$
(b)
$[\mathbf{u}\leftarrow\mathbf{v},\mathbf{u}_{1}\leftarrow\mathbf{v}_{1}^{\prime},\mathbf{u}_{2}^{*}\leftarrow\mathbf{v}_{2}^{*}]\psi$,
for any $\mathbf{u}_{2}^{*}\subseteq\mathbf{u}_{2}$ and $\mathbf{v}_{2}^{*}$
is the original factual value of $\mathbf{u}_{2}^{*}$ (a subvector of
$\mathbf{v}_{2}$).

  3. 3.

AC3: $\nexists\mathbf{u}{\color[rgb]{0,0,1}{}^{\prime}}\subset\mathbf{u}$ that
satisfies AC1 and AC2.

Condition AC1 ensures the suspected actual cause and outcome are factual.
Then, AC2(a) ensures a sufficient counterfactual, that switching the suspect
$\mathbf{u}$ and a circumstance $\mathbf{u}_{1}$ to that counterfactual
assignment guarantees a flipped outcome $\neg\psi$, and AC2(b) states that
switching the circumstance alone does not change the outcome - as long as the
suspect remains the factual value assignments. Finally, AC3 guarantees
minimality of the actual cause. Detailed explanation for this definition can
be found in the original publications [3, 4].

Researchers have applied actual causality and HP model to CPS diagnosis. For
instance, Ibrahim et al. have designed a SAT solver to practically compute
actual causes to explain undesirable CPS behaviors [29, 30]. However, the
solver is restricted to finite $\mathcal{V}(u)$ for variables and their HP
model design only captures discrete events like ”there exists a Byzantine
fault” or ”the system is on autopilot mode”. Unfortunately, this design does
not extend to diagnosis of LECs, which generally have continuous value spaces
for I/O and internal variables, and their internal information flows are not
interpretable.

##  III Problem Formulation

###  III-A System Setting

We formalize the system setting as follows. We have an exact dynamical system
model $\mathcal{S}$, which describes how the components of the agent interact
with each other, as well as how the agent interacts with the environment. This
closed-loop model is assumed to be deterministic. In other words, the agent’s
trajectory only depends on the initial state, the environment dynamics, and
the component designs. These assumptions are often satisfied in model-based
CPS engineering and can be relaxed in future research.

The system aims to satisfy a given specification/property $\varphi$ that takes
Boolean values (true/satisfied/1, false/violated/0), e.g., the linear temporal
logic (LTL) [31] or signal temporal logic (STL) [32] formula. In this paper,
we focus on STL.

We observe a trajectory on a given initial state $s_{0}$, where the system
failed to satisfy $\varphi$ and we suspect one of its components
$\mathcal{C}$, such as a controller, is at fault for the violation. The
component $\mathcal{C}$’s internal design is invisible to us, but we can
observe its input/output (I/O) behavior, which we denote by function
$f:\mathcal{X}\mapsto\mathcal{Y}$, with its domain and codomain being
continuous and bounded metric spaces $(\mathcal{X},d_{\mathcal{X}})$ and
$(\mathcal{Y},d_{\mathcal{Y}})$, respectively. Therefore, we can replace this
behavior by any counterfactual $f^{\prime}:\mathcal{X}\mapsto\mathcal{Y}$. We
denote the counterfactual system that uses $f^{\prime}$ in place of $f$, with
everything else remaining the same, as $\mathcal{S}(f^{\prime})$. We use
$\mathcal{S}(f^{\prime})\models\varphi$ and
$\mathcal{S}(f^{\prime})\not\models\varphi$ to denote that the property will
be satisfied or violated under $\mathcal{S}(f^{\prime})$. With this notation,
the factual outcome is $\mathcal{S}(f)\not\models\varphi$. Next, we formalize
the distance between two choices for the behaviors of component $\mathcal{C}$.

######  Definition 2 (Distance between I/O behaviors).

For any two I/O behaviors $f_{1},f_{2}:\mathcal{X}\mapsto\mathcal{Y}$, we
define distance $||\cdot||_{d_{f}}$ as

$||f_{1}-f_{2}||_{d_{f}}=\max_{x\in\mathcal{X}}||f_{1}(x)-f_{2}(x)||_{d_{\mathcal{Y}}}$
(1)

We then make the following assumptions.

######  Assumption 1.

We can check the outcome of $\varphi$ on the given initial state $s_{0}$ when
substituting different $f^{\prime}$ in place of the component by calling a
simulator,

$\textsc{simulate}_{s_{0}}:(\mathcal{X}\mapsto\mathcal{Y})\mapsto\\{0,1\\},$
(2)

which encodes the knowledge of $\mathcal{S}$ and $\varphi$. For simplicity, we
assume a fixed $s_{0}$ in the remainder of this paper, and the simulator is
denoted as simply simulate.

######  Assumption 2.

The behavior $f$ is Lipschitz-continuous, with an unknown Lipschitz constant.
This is a common property of many types of learning models, such as neural
networks.

######  Assumption 3.

The property outcome is robust against small changes from the factual I/O
behavior $f$, i.e., $\forall f^{\prime}:\mathcal{X}\mapsto\mathcal{Y}$,

$||f-f^{\prime}||_{d_{f}}\leq\epsilon\implies\textsc{simulate}(f)=\textsc{simulate}(f^{\prime}),$
(3)

for some small $\epsilon>0$.

We expect Assumption 3 to hold in most practical cases. Suppose the distance
of the factual $f$ to the decision boundary of the simulator outcome on I/O
behaviors is $\delta$. If $\delta>0$, there exists an arbitrarily small
$\epsilon<\delta$ where the assumption holds. The only case that Assumption 3
does not hold is when $\delta=0$, i.e., the factual behavior is right on the
decision boundary, but this event would usually have a probability measure of
$0$.

###  III-B Problem Statement

In the above setting, we want to identify a subset of I/O behaviors of the
suspected component $\mathcal{C}$ — that is, a set of input-output tuples of
$f$ — that indeed caused the property violation, as well as the counterfactual
outputs on these inputs that can repair the system.

Main Problem. Upon the observation of a property violation on a runtime trace,
$\mathcal{S}(f)\not\models\varphi$, how can we use HP causality to identify a
subset of the suspected component $\mathcal{C}$’s I/O behaviors (modeled as
$f$) that caused this violation?

  1. 1.

Sub-problem 1. Encode the dependency structure of $\mathcal{C}$’s behaviors on
$\varphi$ using an HP model.

  2. 2.

Sub-problem 2. On the encoded HP model, design a search algorithm for an
actual cause, such that, upon success, provide a repair suggestion in form of
a counterfactual $f^{*}:\mathcal{X}\mapsto\mathcal{Y}$ that
$\mathcal{S}(f^{*})\models\varphi$. Upon failure, quantify the confidence that
the property violation is not caused by the I/Os of $\mathcal{C}$.

To solve these problems, we employ the following workflow, which is visualized
in Figure 1: (1) Extract the behaviors of $\mathcal{C}$ as an I/O table. (2)
Encode the dependency of the property outcome on the I/O behaviors with an HP
model (Section IV). (3) Search for a counterfactual model value assignment,
revealing an actual cause and a repair (Section V). (4) Decode the found
assignment as a counterfactual component behavior. (5) Replace $\mathcal{C}$
with an alternative component that performs this counterfactual behavior to
repair the system.

The entire workflow is implemented on an OpenAI Gym example in Section VI.

Figure 1: Workflow of our causal repair approach, which constructs an HP model
to encode the I/O behaviors of a suspected LEC and search for a repair by
causal analysis.

##  IV Halpern-Pearl Model Design

We first describe a naive way of encoding the dependency of $\varphi$ on
component $\mathcal{C}$’s I/O behaviors in Section IV-A. This results in an HP
model with an infinite number of nodes. To alleviate this issue, we describe a
method for constructing an HP model with a finite number of nodes in Section
IV-B.

To fulfill the minimality in AC3 of Definition 1, we will need to distinguish
which counterfactuals are closer to the factual behavior. This requires a
partial order on I/O behaviors that aligns with some partial order on the HP
node values. Unfortunately, the naive finite HP model does not admit fine
enough partial orders over its node values. Thus, we will refine our HP model
to support suitable partial orders and make it amenable to the counterfactual
search.

###  IV-A Infinite HP Model

Intuitively, a naive way to build the HP model to model the dependency of
property outcome on I/Os of $\mathcal{C}$ is to model every input
$x\in\mathcal{X}$ as an endogenous node, and its corresponding output
$y\in\mathcal{Y}$ as its value. We illustrate this model in Figure 2 and refer
to it as the _infinite HP model_.

Figure 2: Infinite HP model.

This model encodes the component $\mathcal{C}$’s implemented I/O behaviors as
an exogenous node (yellow), and the behaviors for each input
$x\in\mathcal{X}$, as well as the property outcome as endogenous nodes. Since
we assume the full knowledge of the $\mathcal{C}$’s I/Os, we can extract the
input-output tuples by tabulating them. Then, based on Assumption 1, the
outcome can be obtained by calling simulate on the I/Os. This design is
illustrated in Figure 2. Each node has its variable name at the top and a
value at the bottom.

Recall that the input space $\mathcal{X}$ is continuous. Therefore, an obvious
drawback of this infinite HP model is that it has infinitely many nodes for
each $x\in\mathcal{X}$, and thus an infinitely large search space for an
actual cause given by Definition 1.

###  IV-B Discretized HP Model

To create a finite HP model, we partition the input space $\mathcal{X}$ and
output space $\mathcal{Y}$ into sufficiently fine-grained, but finitely many
cells $\mathcal{X}=\bigcup_{i=1}^{m}\mathcal{X}_{i}$ and
$\mathcal{Y}=\bigcup_{j=1}^{n}\mathcal{Y}_{j}$, and then consider functions
that map all $x\in\mathcal{X}_{1}$ to the center of one $\mathcal{Y}_{j}$. The
partitioning procedure is detailed in Algorithm 1.

Input : Factual behavior $f:\mathcal{X}\mapsto\mathcal{Y}$, continuous bounded
spaces $\mathcal{X}$ and $\mathcal{Y}$, initial cell width
$\Delta_{x}^{\text{init}}$, $\Delta_{y}^{\text{init}}$

Output : Partitioned cells $\mathcal{X}_{1},\dots,\mathcal{X}_{m}$,
$\mathcal{Y}_{1},\dots,\mathcal{Y}_{n}$, a map
$g:\\{1,\dots,m\\}\mapsto\\{1,\dots,n\\}$ from input cells to output cells

1 $\Delta_{x}\leftarrow\Delta_{x}^{\text{init}}$;

2 $\Delta_{y}\leftarrow\Delta_{y}^{\text{init}}$;

3 do

4 Partition $\mathcal{Y}$ with cell width $\Delta_{y}$;

5 $\Delta_{y}\leftarrow\Delta_{y}/2$;

6 do

7 Partition $\mathcal{X}$ with cell width $\Delta_{x}$;

8 $\Delta_{x}\leftarrow\Delta_{x}/2$;

9

10 while  _$\exists\mathcal{X}_{i}$ that cannot be completely mapped into some
$\mathcal{Y}_{j}$_;

11 $g\leftarrow$ the current mapping between cells;

12 $f_{r}\leftarrow\textsc{recon}(g)$;

13

14while  _$\textsc{simulate}(f_{r})\neq\textsc{simulate}(f)$_ ;

Algorithm 1 Discretize HP Model

Algorithm 1 splits the input space and output space into hypercube cells. The
cell size keeps shrinking until the two conditions, respectively in line 10
and line 13 are met. We define the reconstruction method recon at line 12 as
follows.

$\textsc{recon}(g)=f_{r},\text{ where }\forall
x\in\mathcal{X}_{i},f_{r}(x)=\text{center}(\mathcal{Y}_{j}),j=g(i)$

That is, it reconstructs an I/O behavior from the discrete map $g$, which
approximates the factual I/O behaviors as a mapping among cells, by mapping
every $x$ in an input cell to the center of the corresponding output cell.
This algorithm is guaranteed to terminate by the following theorem.

######  Theorem 1 (Termination of Discretization).

Algorithm 1 is guaranteed to terminate.

For the proof, please refer to Appendix A. Proofs of later theorems can also
be found there.

Upon the termination of Algorithm 1, we obtain a function
$f_{r}=\textsc{recon}(g)$ that represents $f$ by (1) having arbitrarily close
outputs on the same inputs and (2) having the same property outcome. Then, our
causal analysis is established on the family of functions like $f_{r}$, i.e.,

######  Definition 3 (Representative Component Behavior Space).

Given the cell partition $\mathcal{X}=\bigcup_{i}^{m}\mathcal{X}_{i}$ and
$\mathcal{Y}=\bigcup_{j}^{n}\mathcal{Y}_{j}$, we define the _representative
component behavior space_ as a finite subset of
$(\mathcal{X}\mapsto\mathcal{Y})$ as

$\mathcal{F}_{r}=\\{f^{\prime}\in(\mathcal{X}\mapsto\mathcal{Y})\mid\forall
i,\forall x\in\mathcal{X}_{i},\exists
j,f^{\prime}(x)=\operatorname{center}(\mathcal{Y}_{j})\\}$ (4)

Notice that if we wrap up the generation of $g$ from $f$ in Algorithm 1 as a
method discretize, the two methods discretize and recon are inverse to each
other if the behaviors are restricted within $\mathcal{F}_{r}$.

With this shrunken space of behavior choices, we modify our HP model to the
next version, called the _discretized HP model_ , as in Figure 3. Here, we
have one endogenous node per input cell, and its value is the mapped output
cell. This HP model has $m+2$ nodes and can express every behavior choice in
$\mathcal{F}_{r}$. Notice that in place of tabulate, we now have discretize
since the nodes are now representing the discrete map $g$, and the simulate
method requires recon on the discrete map first.

Figure 3: Discretized HP model. We illustrate this with an example value
assignment of nodes.

Generally, we want to pick a counterfactual for $\mathcal{C}$ that repairs the
system $\mathcal{S}$, i.e., flipping the outcome from 0 to 1, that is closest
to the factual. This is motivated by the minimality of actual causes in AC3 in
Definition 1. Therefore, we first define a partial order on I/O behaviors.

######  Definition 4 (Partial Order on Behaviors).

For a factual behavior choice $f:\mathcal{X}\mapsto\mathcal{Y}$, we define a
_partial behavior order_ $\preccurlyeq_{f}$ on
$(\mathcal{X}\mapsto\mathcal{Y})$ as

$\begin{split}f_{1}\preccurlyeq_{f}f_{2}&\Longleftrightarrow\forall
x\in\mathcal{X},\forall j\in\\{1,\dots,\operatorname{dim}(\mathcal{Y})\\},\\\
&(f(x)[j]\leq f_{1}(x)[j]\leq f_{2}(x)[j])\\\ \lor&(f(x)[j]\geq
f_{1}(x)[j]\geq f_{2}(x)[j])\end{split}$ (5)

where $[j]$ denotes the $j$-th dimension of a vector.

In plain words, $f_{1}\preccurlyeq_{f}f_{2}$ iff $f_{1}$ has closer outputs to
the factual $f$ than $f_{2}$ does on all dimensions and on all inputs.

However, one drawback of this discrete HP model and the partial order from
Definition 4 is that we cannot tell the difference between two counterfactual
I/O behaviors in terms of the sets of “disagreeing” nodes compared to the
factual $\mathbf{v}=\textsc{discretize}(f_{r})$ , which is something we need
for reasoning about (AC3) in Definition 1. For example, in Figure 3, assume
that two counterfactuals $f_{1}$ and $f_{2}$ change the mapping on
$\mathcal{X}_{1}$ from the factual $\mathcal{Y}_{2}$ (as in $f_{r}$) to
$\mathcal{Y}_{3}$ and $\mathcal{Y}_{4}$, respectively. With this mapping
change, both flip the outcome to 1. In the discretized HP model, the number of
“disagreeing” nodes under both assignments is 1, but $\mathcal{Y}_{4}$ may be
further from the factual $\mathcal{Y}_{2}$ than $\mathcal{Y}_{3}$. In plain
words, we cannot answer “which one is more different from $f_{r}$: $f_{1}$ or
$f_{2}$?” by simply comparing the two sets of nodes they modify from the
factual value assignment.

###  IV-C Propositional HP Model

We now present a partial order, $\preccurlyeq_{\mathbf{v}}$, over node value
differences and define how $\preccurlyeq_{\mathbf{v}}$ relates to partial
order $\preccurlyeq_{f}$ from Definition 4. We then present an HP model
construction technique, starting from the discretized HP model, that preserves
$\preccurlyeq_{\mathbf{v}}$ and $\preccurlyeq_{f}$.

For the new partial order, if the subset of nodes that differ between
assignments $\mathbf{v}_{2}$ and factual $\mathbf{v}$ contains the subset of
nodes that differ between $\mathbf{v_{1}}$ and $\mathbf{v}$, we say that
$\mathbf{v}_{1}$ is closer to $\mathbf{v}$ than $\mathbf{v}_{2}$ does. We
formulate this as another partial order, on value assignments of a subset of
endogenous nodes $\mathcal{U}_{\text{io}}\subseteq\mathcal{U}_{\text{endo}}$.
The value space of $\mathcal{U}_{\text{io}}$ is denoted as
$\mathcal{V}_{\text{io}}$.

######  Definition 5 (Partial Order on HP Node Values).

On a set of HP model nodes $\mathcal{U}$ and its value space $\mathcal{V}$,
given the factual value assignment $\mathbf{v}\in\mathcal{V}$, we can define a
partial order $\preccurlyeq_{\mathbf{v}}$ that

$\mathbf{v}_{1}\preccurlyeq_{\mathbf{v}}\mathbf{v}_{2}\Longleftrightarrow\operatorname{diff}_{\mathcal{U}}(\mathbf{v},\mathbf{v}_{1})\subseteq\operatorname{diff}_{\mathcal{U}}(\mathbf{v},\mathbf{v}_{2}),$
(6)

where $\operatorname{diff}_{\mathcal{U}}(\cdot,\cdot)$ denotes the subset of
nodes in $\mathcal{U}$ that has different values by two assignments. Equality
$=_{\mathbf{v}}$ holds iff
$\operatorname{diff}_{\mathcal{U}}(\mathbf{v},\mathbf{v}_{1})=\operatorname{diff}_{\mathcal{U}}(\mathbf{v},\mathbf{v}_{2})$.

With this partial order, we can determine if an assignment $\mathbf{v}_{1}$
differs more from the factual $\mathbf{v}$ than assignment $\mathbf{v}_{2}$
does.

Consequently, a larger change from the factual I/O behavior needs to be
reflected in a larger set of differing nodes, i.e., the partial order
$\preccurlyeq_{\mathbf{v}}$ needs to preserve the partial order
$\preccurlyeq_{f}$. We define partial order preservation as follows:

######  Definition 6 (Partial Order Preservation).

Consider a representative behavior space
$\mathcal{F}_{r}\subset(\mathcal{X}\mapsto\mathcal{Y})$ with a factual
representative component $f_{r}$ and the induced partial behavior order
$\preccurlyeq_{f_{r}}$. Consider also a subset of endogenous nodes,
$\mathcal{U}_{\text{io}}\subseteq\mathcal{U}_{\text{endo}}$, with the value
space of this subset $\mathcal{V}_{\text{io}}$. An encoding of a
representative I/O behavior onto the nodes $\mathcal{U}_{\text{io}}$,
$\textsc{encode}:\mathcal{F}_{r}\mapsto\mathcal{V}_{\text{io}}$, _preserves
partial order_ iff $\textsc{encode}(f_{r})=\mathbf{v}$, and

$\forall
f_{1},f_{2}\in\mathcal{F}_{r},f_{1}\preccurlyeq_{f_{r}}f_{2}\Longleftrightarrow\textsc{encode}(f_{1})\preccurlyeq_{\mathbf{v}}\textsc{encode}(f_{2})$
(7)

Finally, we define a finer-grained HP model that preserves these two partial
orders based on Definition 6.

First, we first define an output bin of $\mathcal{Y}$. As illustrated in
Figure 4, if the output space is partitioned into $4\times 4$ hypercube cells,
we can index the cells from $\mathcal{Y}_{11}$ to $\mathcal{Y}_{44}$. We
therefore have 4 bins along dimension 1 and 4 bins along dimension 2.
Formally, if we index the $\text{dim}(\mathcal{Y})=d$-dimensional output cells
as $\mathcal{Y}_{k_{1}k_{2}\dots k_{d}}$, we have

######  Definition 7 (Output Bins).

The $k$-th output bin along dimension $j$ is
$bin_{\mathcal{Y}}(j,k)=\bigcup\\{\mathcal{Y}_{k_{1}k_{2}\dots k_{d}}\mid
k_{j}=k\\}$

We denote the total number of bins along the $j$-th dimension as $n_{j}$, and
therefore $k_{j}\in\\{1,\dots,n_{j}\\}$. Next, we can define the lower bound
of $bin_{\mathcal{Y}}(j,k)$ in $j$-th dimension as
$lo(bin_{\mathcal{Y}}(j,k))$. Notice that along a fixed dimension $j$, we have
a total order of lower bounds of bins based on $k$.

We are now ready to construct the final, _propositional HP model_ , created by
Algorithm 2. In contrast to the discretized HP model, now the endogenous nodes
that encode the I/O behaviors ($\mathcal{U}_{\text{io}}$ in the algorithm)
take propositional values. Node $u_{ijk}$ encodes whether for input
$x\in\mathcal{X}_{i}$ the output cell for the $j$-th dimension is in at least
the $k$-th output bin. In the example illustrated in Figure 5, we have
1-dimensional $\mathcal{X}$, with $\mathcal{X}_{1}=[0,1]$ and
$\mathcal{X}_{2}=[1,2]$. We also have 2-dimensional
$\mathcal{Y}=([5,6]\cup[6,7])\times([10,11]\cup[11,12])$. The example value
shows that the representative function $f_{r}$ of $f$ maps all $x\in[0,1]$ to
the center of $[6,7]\times[10,11]$, and all $x\in[1,2]$ to
$[5,6]\times[10,11]$. Under this propositional HP model, the encoding
$\textsc{encode}:\mathcal{V}_{\text{io}}\mapsto\mathcal{F}_{r}$ is simply
evaluating the node propositions based on the component behavior, and decode
is first identifying the discretized mapping between cells, and then calling
recon. Notice that encode and decode are inverse to each other if we restrict
the behavior choices in $\mathcal{F}_{r}$.

Figure 4: An example indexing of cells and output bins when both $\mathcal{X}$
and $\mathcal{Y}$ are 2-dimensional.

Input : Input space $\mathcal{X}$ and output space $\mathcal{Y}$ of component
behaviors and their indexed cells, simulator simulate, encoder encode, decoder
decode

Output : HP model
$\mathcal{M}=(\mathcal{U}_{\text{endo}},\mathcal{U}_{\text{exo}},\mathcal{V},\mathcal{E})$

1 $\mathcal{U}_{\text{exo}}\leftarrow\\{u_{\text{comp}}$};

2 $\mathcal{U}_{\text{endo}}\leftarrow\\{u_{\varphi}\\}$;

3 $\mathcal{V}(u_{\text{comp}})\leftarrow(\mathcal{X}\mapsto\mathcal{Y})$;

4 $\mathcal{U}_{\text{io}}\leftarrow\emptyset$;

5
$\mathcal{E}\leftarrow\\{\textsc{simulate}\circ\textsc{decode},\textsc{encode}\\}$;

6 for  _$i=1,\dots,m$_ do

7 for  _$j=1,\dots,d$_ do

8 for  _$k=1,\dots,n_{j}$_ do

9 $u_{ijk}\leftarrow(x\in\mathcal{X}_{i}\implies f(x)[k]\geq
lo(bin_{\mathcal{Y}}(j,k))$;

10 $\mathcal{V}(u_{ijk})\leftarrow\\{0,1\\}$ ;

11
$\mathcal{U}_{\text{io}}\leftarrow\mathcal{U}_{\text{io}}\cup\\{u_{ijk}\\}$;

12

13 end for

14

15 end for

16

17 end for

18$\mathcal{U}_{\text{endo}}\leftarrow\mathcal{U}_{\text{io}}\cup\mathcal{U}_{\text{endo}}$;

Algorithm 2 Construct Propositional HP Model Figure 5: An example of the
propositional HP model. We use red and green colors to show the nodes with 0
and 1 values.

With $m$ input cells, $d$ output dimensions and $n_{j}$ bins along a dimension
$j$, this HP model has $m\sum_{j=1}^{d}n_{j}+2$ nodes, with
$|\mathcal{U}_{\text{io}}|=m\sum_{j=1}^{d}n_{j}$. However, the total number of
value assignments on $\mathcal{U}_{\text{io}}$ is not
$2^{|\mathcal{U}_{\text{io}}|}$, because some value assignments are not
allowed. For example, we cannot let $f(x)\geq 5$ be false and $f(x)\geq 6$ be
true. The total number of valid assignments is
$(\prod_{j=1}^{d}n_{j})^{m}=n^{m}$, the same as the discretized HP model.

######  Theorem 2 (Propositional Encoding Preserves Partial Order).

The encoding $\textsc{encode}:\mathcal{F}_{r}\mapsto\mathcal{V}_{\text{io}}$
specified by evaluating the proposition of every
$u_{ijk}\in\mathcal{U}_{\text{io}}$, as constructed in Algorithm 2, preserves
the partial order as defined in Definition 6.

##  V The Causal Repair Algorithm

###  V-A Satisfactory Counterfactual Search

Based on the constructed propositional HP model, we look for an actual cause.
The idea is to first search for a counterfactual
$\mathbf{v}_{\text{io}}^{\prime}$ on $\mathcal{U}_{\text{io}}$ that leads to a
satisfactory outcome, i.e., $\mathcal{S}(f_{r}^{\prime})\models\varphi$ and
$f_{r}^{\prime}=\textsc{decode}(\mathbf{v}_{\text{io}}^{\prime})$. The subset
of nodes with different value assignments between the factual
$\mathbf{v}_{\text{io}}$ and counterfactual $\mathbf{v}_{\text{io}}^{\prime}$
is not necessarily an actual cause yet (it may not be minimal as per condition
AC3), and in the next subsection, we will look for a different
$\mathbf{v}_{\text{io}}^{*}$ that fulfills the actual cause conditions from
this $\mathbf{v}_{\text{io}}^{\prime}$. Here, we focus on the search for
$\mathbf{v}_{\text{io}}^{\prime}$.

Since a brute-force search on all $n^{m}$ possible value assignments takes
exponential time, we leverage random sampling as follows in Algorithm 3.

Input : HP model $\mathcal{M}$, probability threshold $p$, significance level
$\alpha\in[0,1]$, simulator simulate, decoder decode, probability
distributions $\mathcal{D}_{1}$, $\mathcal{D}_{2}$

Output : Either a counterfactual value assignment
$\mathbf{v}_{\text{io}}^{\prime}$ on $\mathcal{U}_{\text{io}}$ with
$\textsc{simulate}(\textsc{decode}(\mathbf{v}_{\text{io}}^{\prime}))=1$, or a
statement
$\Pr_{\mathcal{D}_{2}}[\Pr_{\mathcal{D}_{1}}[\textsc{simulate}(\textsc{decode}(\mathbf{v}_{\text{io}}^{\prime}))=1]\leq
p]\geq 1-\alpha$

1 $N\leftarrow\lceil(1/p-1)Q(1-\alpha/2)^{2}\rceil$;

2 for  _$1\dots N$_ do

3 $\mathbf{v}_{\text{io}}^{\prime}\leftarrow$ uniform sampling from all
$n^{m}$ settings;

4
$\varphi\leftarrow\textsc{simulate}(\textsc{decode}(\mathbf{v}_{\text{io}}^{\prime}))$;

5 if  _$\varphi$_ then

6 Return $\mathbf{v}_{\text{io}}^{\prime}$ ;

7

8 end if

9

10 end for

11Return
$\Pr_{\mathcal{D}_{2}}[\Pr_{\mathcal{D}_{1}}[\textsc{simulate}(\textsc{decode}(\mathbf{v}_{\text{io}}^{\prime}))=1]\leq
p]\geq 1-\alpha$;

Algorithm 3 Counterfactual Random Sampling

The distribution $\mathcal{D}_{1}$ is a distribution on different value
assignments in the value space $\mathcal{V}_{\text{io}}$. For example, it can
be a uniform distribution on the propositional node values. Distribution
$\mathcal{D}_{2}$ is the induced distribution on a sampled value assignment’s
success rate, i.e.,
$p=\Pr_{\mathcal{D}_{1}}[\textsc{simulate}(\textsc{decode}{(\mathbf{v}_{\text{io}}^{\prime})}=1]$
, which is treated as another random variable. In line 1, the function
$Q(\cdot)$ means the quantile on standard normal distribution. Starting from
line 2 to 8, we uniformly sample different value assignments
$\mathbf{v}_{\text{io}}^{\prime}$ until we find one
$\mathbf{v}_{\text{io}}^{\prime}$ that produces satisfactory $\varphi$ by
simulate or we reach a maximum number of samples. If we fail to find a
satisfactory value assignment, we report that the probability of finding a
satisfactory $\mathbf{v}_{\text{io}}^{\prime}$ with uniform sampling, i.e.,
the portion of satisfactory $\mathbf{v}_{\text{io}}^{\prime}$ in the entire
search space, is at most $p$ with confidence $1-\alpha$ at line 9. We show our
probabilistic failure statement at line 9 holds with the following theorem.

######  Theorem 3 (Probabilistic Guarantee on Search Failure).

Given an HP model constructed by Algorithm 2, a probability threshold
$p\in[0,1]$ and a confidence $1-\alpha\in[0,1]$, the final statement at line 9
of Algorithm 3 holds.

Upon $N$ consecutive failures, this search algorithm states that with some
confidence a counterfactual assignment that flips the outcome is unlikely to
be found. This suggests that the actual cause for the run-time property
violation lies elsewhere: possibly the environment, the behaviors of another
component, or the suspected component together with another component — but
not the I/O behaviors of this component alone. This confidence is based on a
substantial number of samples without finding a successful counterfactual. For
example, for $p=0.001$ and $\alpha=0.05$, one would need to uniformly sample
at least $N=3838$ failed counterfactuals in a row.

###  V-B Node Value Interpolation for Actual Cause

Suppose we have successfully obtained a satisfactory
$\mathbf{v}_{\text{io}}^{\prime}$ on $\mathcal{U}_{\text{io}}$, from Algorithm
3. The final step is to find a counterfactual $\mathbf{v}_{\text{io}}^{*}$
that can flip the outcome and is as close to the factual
$\mathbf{v}_{\text{io}}$ as possible, starting from
$\mathbf{v}_{\text{io}}^{\prime}$. This step is required to satisfy the
minimiality condition (AC3) of Definition 1.

We therefore perform a deterministic interpolation between
$\mathbf{v}_{\text{io}}^{\prime}$ and $\mathbf{v}_{\text{io}}$ in Algorithm 4,
which starts from the found satisfactory counterfactual assignment
$\mathbf{v}_{\text{io}}^{\prime}$. In this algorithm, the counterfactual value
assignment steps towards the factual $\mathbf{v}_{\text{io}}$ by flipping the
differing value assignments one-by-one. This procedure continues until there
are no nodes it can flip while still satisfying $\varphi$. Because the total
number of nodes in $\mathcal{U}_{\text{io}}$ is $m\sum_{j=1}^{d}n_{j}$,
Algorithm 4 is guaranteed to output an actual cause with a satisfactory
counterfactual within $O(m\sum_{j=1}^{d}n_{j})$ time. The complexity is linear
in terms of $m$ and $n$, allowing for efficient computation.

We visualize the process of searching for $\mathbf{v}_{\text{io}}^{*}$ as
Figure 6. Without loss of generality, consider only the top left input cell.
The factual behavior $f=\textsc{decode}(\mathbf{v}_{\text{io}})$ maps inputs
in this cell to some output cell and this behavior produces a property
violation. The sampling by Algorithm 3 finds a satisfactory counterfactual
$f^{\prime}=\textsc{decode}(\mathbf{v}_{\text{io}}^{\prime})$, which is an
alternative mapping. Next, the interpolation Algorithm 4 flips the nodes
disagreed by $\mathbf{v}_{\text{io}}$ and $\mathbf{v}_{\text{io}}^{\prime}$
one-by-one towards $\mathbf{v}_{\text{io}}$. This flipping is equivalent to
stepping through the output cells one by one. Eventually, the algorithm
reaches a node assignment where any steps towards $f$ result in $\varphi$
becoming violated, at which point it returns the current node assignment.
Depending on the dimensions, we can take different paths in stepping, i.e. we
can end up either in cell (1) or cell (2), from interpolating in the vertical
or horizontal dimension first, respectively. Mapping the input cell to either
of these two cells represents a valid I/O behavior of $f^{*}$.

Input : HP model $\mathcal{M}$, factual $\mathbf{v}_{\text{io}}$, satisfactory
counterfactual $\mathbf{v}_{\text{io}}^{\prime}$, simulate function simulate

Output : Satisfactory counterfactual $\mathbf{v}_{\text{io}}^{*}$ such that
its difference from $\mathbf{v}_{\text{io}}$ is an actual cause as per
$\mathcal{M}$

1 $\mathbf{v}_{\text{io}}^{*}\leftarrow\mathbf{v}_{\text{io}}^{\prime}$;

2
$\operatorname{diff}_{\mathcal{U}_{\text{io}}}^{*}\leftarrow\operatorname{diff}_{\mathcal{U}_{\text{io}}}(\mathbf{v}_{\text{io}},\mathbf{v}_{\text{io}}^{\prime})$;

3 for  _$i=1,\dots,m$_ do

4 for  _$j=1,\dots,d$_ do

5 for  _$k=1,\dots,n_{d}$_ do

6 if  _$u_{ijk}\notin\operatorname{diff}_{\mathcal{U}_{\text{io}}}^{*}$_ then

7 Continue;

8

9 end if

10 $\mathbf{v}_{\text{temp}}\leftarrow$ $\mathbf{v}_{\text{io}}^{*}$ with
value assignment on $u_{ijk}$ the negation as in $\mathbf{v}_{\text{io}}^{*}$;

11 if  _$\textsc{simulate}(\textsc{decode}(\mathbf{v}_{\text{temp}}))=$ 1_
then

12 $\mathbf{v}_{\text{io}}^{*}\leftarrow\mathbf{v}_{\text{temp}}$;

13
$\operatorname{diff}_{\mathcal{U}_{\text{io}}}^{*}\leftarrow\operatorname{diff}_{\mathcal{U}_{\text{io}}}^{*}\setminus\\{u_{ijk}\\}$

14 end if

15

16 end for

17

18 end for

19

20 end for

Algorithm 4 Actual Cause Search by Incremental Interpolation

Notice that Algorithm 4 incrementally steps towards the factual
$\mathbf{v}_{\text{io}}$ from $\mathbf{v}_{\text{io}}^{\prime}$ by flipping
nodes one-by-one. A more efficient variant is to do a _binary search_ between
these two value assignments. We denote these two approaches as incremental
interpolation and binary search interpolation, respectively. Both are
evaluated in Section VI.

Next, we show that the output is indeed an actual cause based on HP model
$\mathcal{M}$ in Theorem 4.

######  Theorem 4 (Output is Actual Cause).

Let HP model $\mathcal{M}$ constructed in Algorithm 2 be given with factual
node value assignment $\mathbf{v}_{\text{io}}$. Let
$\mathbf{v}_{\text{io}}^{\prime}$ be a counterfactual node value assignment
from Algorithm 3. The node values on $\mathcal{U}_{\text{io}}$ where
assignments $\mathbf{v}_{\text{io}}$ and $\mathbf{v}_{\text{io}}^{*}$ disagree
in Algorithm 4 are an _actual cause_ of $\mathcal{S}(f)\not\models\varphi$ as
per $\mathcal{M}$ constructed in Algorithm 2.

Figure 6: An example of finding actual cause by interpolation, with two-
dimensional $\mathcal{X}$ and $\mathcal{Y}$. Red cells in $\mathcal{Y}$ denote
property violation, green denotes satisfaction, and blue denotes cells where
satisfaction is not relevant to the example.

##  VI Experimental Evaluation

###  VI-A Setup

Figure 7: Control functions, approximated by mapping between cells, visualized
as heatmaps on the position-velocity space. The functions from left to right
are the factual $f$, the searched counterfactual $f_{1}^{\prime}$ and
interpolated counterfactual $f^{*}_{1}$ between the first two.

We test our diagnosis on the mountain car from the OpenAI Gym [13], which has
become a common benchmark for learning-enabled CPS [33, 34]. The car starts in
the valley between two mountains with the task is to drive it to the top of
the right mountain within a deadline. The system has two one-dimensional state
variables: position and velocity, and a one-dimensional control signal, as
follows.

$\begin{split}&pos(t+1)=pos(t)+vel(t)\\\
&vel(t+1)=vel(t)+0.0015ctrl(t)-z\operatorname{cos}{(3pos(t))}\\\
&ctrl(t)=f(pos(t),vel(t))\end{split}$ (8)

where $z=0.0025$ is the steepness of the hill. The variables are bounded, such
that $pos(t)\in[-1.2,0.6]$, $vel(t)\in[-0.07,0.07]$ and $ctrl(t)\in[-1,1]$.
The initial condition is $(pos(0),vel(0))=(-0.5,0)$, i.e., staying still at
the bottom of the valley. The controller is a learning-based function $f$, for
which we use a pre-trained deep neural network. The run-time property can be
specified as an STL [32] formula

$\varphi=F_{t\leq 110}(pos(t)\geq 0.45),$ (9)

i.e., reaching $pos(t)\geq 0.45$ before the deadline of $t=110$. Due to its
limited power, the car cannot reach its goal directly, and the challenge is to
first climb the left mountain to gain enough momentum.

As the controller function $f$, we use a pre-trained rectangular deep neural
network with shape $8\times 16$ with sigmoid activations. The run-time
property $\varphi$ is violated on initial state $(pos(0),vel(0))=(-0.5,0)$ ,
and we suspect it is caused by this learned controller. Therefore, we run our
diagnosis step-by-step to find an actual cause and observe its effect on
system repair. All experiments are run on a single core of Intel Xeon Gold
6148 CPU @ 2.40GHz.

###  VI-B Results

Based on the setup, the input space and output space of the suspected
controller are $\mathcal{X}=[-0.6,1.2]\times[-0.07,0.07]$ and
$\mathcal{Y}=[-1,1]$. We first run Algorithm 1 and find that the cell width is
of size $0.1\times 0.01$ on $\mathcal{X}$ and $0.1$ on $\mathcal{Y}$. In other
words, position is split into $(1.2-(-0.6))/0.1=18$ equally sized intervals,
velocity into $(0.07-(-0.07))/0.01=14$ intervals and control into
$(1-(-1))/0.1=20$ intervals. Therefore, there are $18\times 14=252$ input
cells and $20$ output cells.

By using Algorithm 2, we constructed an HP model, with every
$u_{ik}\in\mathcal{U}_{\text{io}}$ representing a Boolean statement “on this
input cell $\mathcal{X}_{i}$ of (position, velocity), the output control is at
least in output bin $bin_{\mathcal{Y}}(1,k)$”. We do not have $j$ here because
the output space is one-dimensional. Then, using Algorithm 3, we found three
counterfactuals $\mathbf{v}_{\text{io}}^{1}$, $\mathbf{v}_{\text{io}}^{2}$,
and $\mathbf{v}_{\text{io}}^{3}$ that lead to the satisfaction of $\varphi$.
We then applied Algorithm 4 to find the minimal modification needed from
factual $\mathbf{v}_{\text{io}}$ to these 3 counterfactuals, with the
interpolated counterfactuals denoted as $\mathbf{v}_{\text{io}}^{1*}$,
$\mathbf{v}_{\text{io}}^{2*}$ and $\mathbf{v}_{\text{io}}^{3*}$, respectively.

Figure 7 shows the control functions $f$,
$f^{\prime}_{1}=\textsc{decode}(\mathbf{v}_{\text{io}}^{1})$, and
$f^{*}_{1}=\textsc{decode}(\mathbf{v}_{\text{io}}^{1*})$ as heatmaps. There
are 153 (out of 252) input cells that map to a different output cell between
$f$ and $f^{*}_{1}$. Consequently, our diagnosis pipeline concludes that the
actual cause is _“the control signals on these 153 input cells given by the
factual controller $f$”_, with a corresponding repair _“had these 153 cells
been mapped by $f^{*}_{1}$ instead of $f$, the system would have satisfied
$\varphi$”_.

After we found a counterfactual behavior, we used both the incremental and
binary search interpolations. The computation time is listed in Table I. We
can see the time overhead is predominantly from running the simulator, and
that interpolation with binary search does reduce this overhead.

Interpolation | Total time (s) | Simulator time (s) | Stepping time (s) | # of operations
---|---|---|---|---
Incremental | 9148.77 | 9148.76 | 0.003 | 1231
Binary | 4965.03 | 4965.02 | 0.001 | 880
TABLE I: Computation time of interpolation from counterfactual
$f_{1}^{\prime}$ to factual $f$. An “operation” is a combination of executing
simulate and stepping the counterfactual towards the factual behavior.

To validate that these modifications indeed repair the system, we replace the
controller with interpolated control functions $f^{*}_{1}$, $f^{*}_{2}$ and
$f^{*}_{3}$. The re-runs give satisfactory results as shown in Figure 8,
showing the utility of our causal diagnosis.

Figure 8: Mountain car traces: one before repair and three causally repaired.
The red dashed line means $pos=0.45$, i.e., the target. The car stops
immediately after target is reached.

##  VII Discussion and Conclusion

Our diagnosis and repair pipeline has several advantages. First, we are able
to leverage domain knowledge of the system, such as how to simulate the
outcome of the given STL property, and well-defined actual causality models to
reason about the causes. Second, the actual cause we identify reflects a
minimal subset of the component’s I/O behaviors that cause the problem under
our defined partial order. Therefore, we avoid changing the complex and
possibly unavailable internal structure of components and blaming unrelated,
“innocent” behaviors. Finally, as the experiments demonstrate, our approach
produces practically useful repair suggestions.

Limitations do exist in our current pipeline, so future research is required.
First, the HP model construction and component behavior encoding in Section
IV-C is only one way of preserving the partial order on behaviors. There may
be more effective ways with fewer HP model nodes and therefore smaller memory
and time overhead. Second, the search algorithm in Section V-A draws uniformly
random samples from an exponentially large search space, and even if we are
confident to report the satisfactory portion is relatively small and hence the
search fails, that portion may still be large in an absolute sense. For
example, a fraction of $10^{-6}$ on a space size of $10^{20}$ is still very
large. Therefore, there is potential for making the search algorithm smarter,
such as leveraging prior knowledge of the STL property. Moreover, we perform a
repair on only one initial condition, which may violate the property on the
others.

To summarize, this paper uses a novel construct of the Halpern-Pearl model to
search for an actual cause for a run-time CPS property violation and provide a
repair, without changing complex internal structures of LECs. Our diagnostic
pipeline consists of two steps: (1) construct the propositional HP model and
(2) search for different variable assignments on the model until we reach an
actual cause. Our causal pipeline provides repair guarantees and can generate
useful fixes to deep neural network controllers.

## Acknowledgement

This work was supported in part by ARO W911NF-20-1-0080 and AFRL and DARPA
FA8750-18-C-0090. Any opinions, findings and conclusions or recommendations
expressed in this material are those of the authors and do not necessarily
reflect the views of the Air Force Research Laboratory (AFRL), the Army
Research Office (ARO), the Defense Advanced Research Projects Agency (DARPA),
or the Department of Defense, or the United States Government.

## References

  * [1] S. Wang, A. Ayoub, B. Kim, G. Gössler, O. Sokolsky, and I. Lee, “A causality analysis framework for component-based real-time systems,” in _International Conference on Runtime Verification_. Springer, 2013, pp. 285–303.
  * [2] J. Schumann, P. Moosbrugger, and K. Y. Rozier, “R2u2: monitoring and diagnosis of security threats for unmanned aerial systems,” in _Runtime Verification_. Springer, 2015, pp. 233–249.
  * [3] J. Y. Halpern and J. Pearl, “Causes and explanations: A structural-model approach. part i: Causes,” _The British journal for the philosophy of science_ , 2005.
  * [4] J. Y. Halpern, _Actual causality_. MiT Press, 2016.
  * [5] H. Chockler and J. Y. Halpern, “Responsibility and blame: A structural-model approach,” _Journal of Artificial Intelligence Research_ , vol. 22, pp. 93–115, 2004.
  * [6] D. Held, S. Thrun, and S. Savarese, “Learning to track at 100 fps with deep regression networks,” in _European conference on computer vision_. Springer, 2016, pp. 749–765.
  * [7] C. Hartsell, N. Mahadevan, S. Ramakrishna, A. Dubey, T. Bapty, T. Johnson, X. Koutsoukos, J. Sztipanovits, and G. Karsai, “Model-based design for cps with learning-enabled components,” in _Proceedings of the Workshop on Design Automation for CPS and IoT_ , 2019, pp. 1–9.
  * [8] C. E. Tuncali, J. Kapinski, H. Ito, and J. V. Deshmukh, “Reasoning about safety of learning-enabled components in autonomous cyber-physical systems,” in _Proceedings of the 55th Annual Design Automation Conference_ , 2018, pp. 1–6.
  * [9] T. Chakraborti, S. Sreedharan, Y. Zhang, and S. Kambhampati, “Plan explanations as model reconciliation: Moving beyond explanation as soliloquy,” _arXiv preprint arXiv:1701.08317_ , 2017.
  * [10] B. Krarup, M. Cashmore, D. Magazzeni, and T. Miller, “Model-based contrastive explanations for explainable planning,” 2019.
  * [11] V. Raman and H. Kress-Gazit, “Explaining impossible high-level robot behaviors,” _IEEE Transactions on Robotics_ , vol. 29, no. 1, 2012.
  * [12] “Towards minimal explanations of unsynthesizability for high-level robot behaviors,” in _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2013, pp. 757–762.
  * [13] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, “Openai gym,” _arXiv preprint arXiv:1606.01540_ , 2016.
  * [14] X. Qi, Y. Luo, G. Wu, K. Boriboonsomsin, and M. Barth, “Deep reinforcement learning enabled self-learning control for energy efficient driving,” _Transportation Research Part C: Emerging Technologies_ , vol. 99, pp. 67–81, 2019.
  * [15] D. Cao, J. Zhao, W. Hu, N. Yu, F. Ding, Q. Huang, and Z. Chen, “Deep reinforcement learning enabled physical-model-free two-timescale voltage control method for active distribution systems,” _IEEE Transactions on Smart Grid_ , vol. 13, no. 1, pp. 149–165, 2021.
  * [16] R. McConnell and R. Rajsuman, “Test and repair of large embedded drams. i,” in _Proceedings International Test Conference 2001 (Cat. No. 01CH37260)_. IEEE, 2001, pp. 163–172.
  * [17] Y. Zorian, “Embedded memory test and repair: Infrastructure ip for soc yield,” in _Proceedings. International Test Conference_. IEEE, 2002, pp. 340–349.
  * [18] D. Cazes and M. Kalech, “Model-based diagnosis with uncertain observations,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 03, 2020, pp. 2766–2773.
  * [19] E. Schulte, J. DiLorenzo, W. Weimer, and S. Forrest, “Automated repair of binary and assembly programs for cooperating embedded devices,” _ACM SIGARCH Computer Architecture News_ , vol. 41, no. 1, pp. 317–328, 2013.
  * [20] U. S. Cruz, J. Ferlez, and Y. Shoukry, “Safe-by-repair: a convex optimization approach for repairing unsafe two-level lattice neural network controllers,” _arXiv preprint arXiv:2104.02788_ , 2021.
  * [21] R. Peddi and N. Bezzo, “Interpretable run-time prediction and planning in co-robotic environments,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 2504–2510.
  * [22] M. J. Islam, R. Pan, G. Nguyen, and H. Rajan, “Repairing deep neural networks: Fix patterns and challenges,” in _Proceedings of the ACM/IEEE 42nd International Conference on Software Engineering_ , 2020, pp. 1135–1146.
  * [23] K. Majd, S. Zhou, H. B. Amor, G. Fainekos, and S. Sankaranarayanan, “Local repair of neural networks using optimization,” _arXiv preprint arXiv:2109.14041_ , 2021.
  * [24] J. Sohn, S. Kang, and S. Yoo, “Search based repair of deep neural networks,” _arXiv preprint arXiv:1912.12463_ , 2019.
  * [25] M. Usman, D. Gopinath, Y. Sun, Y. Noller, and C. S. Păsăreanu, “Nn repair: constraint-based repair of neural network classifiers,” in _Computer Aided Verification: 33rd International Conference, CAV 2021, Virtual Event, July 20–23, 2021, Proceedings, Part I 33_. Springer, 2021, pp. 3–25.
  * [26] X. Lin, H. Zhu, R. Samanta, and S. Jagannathan, “Art: abstraction refinement-guided training for provably correct neural networks,” in _Formal Methods in Computer-aided Design_ , vol. 1. TU Wien Academic Press, 2020, pp. 148–157.
  * [27] D. Cohen and O. Strichman, “Automated repair of neural networks,” _arXiv preprint arXiv:2207.08157_ , 2022.
  * [28] F. Fu, Z. Wang, J. Fan, Y. Wang, C. Huang, X. Chen, Q. Zhu, and W. Li, “Reglo: Provable neural network repair for global robustness properties,” in _Workshop on Trustworthy and Socially Responsible Machine Learning, NeurIPS 2022_.
  * [29] A. Ibrahim, S. Rehwald, and A. Pretschner, “Efficient checking of actual causality with sat solving,” _Engineering Secure and Dependable Software Systems_ , vol. 53, p. 241, 2019.
  * [30] A. Ibrahim, S. Kacianka, A. Pretschner, C. Hartsell, and G. Karsai, “Practical causal models for cyber-physical systems,” in _NASA Formal Methods Symposium_. Springer, 2019, pp. 211–227.
  * [31] A. Pnueli, “The temporal logic of programs,” in _18th Annual Symposium on Foundations of Computer Science (sfcs 1977)_. ieee, 1977, pp. 46–57.
  * [32] O. Maler and D. Nickovic, “Monitoring temporal properties of continuous signals,” in _Formal Techniques, Modelling and Analysis of Timed and Fault-Tolerant Systems_. Springer, 2004, pp. 152–166.
  * [33] R. Ivanov, T. Carpenter, J. Weimer, R. Alur, G. Pappas, and I. Lee, “Verisig 2.0: Verification of Neural Network Controllers Using Taylor Model Preconditioning,” in _Computer Aided Verification_. Cham: Springer International Publishing, 2021, pp. 249–262.
  * [34] I. Ruchkin, M. Cleaveland, R. Ivanov, P. Lu, T. Carpenter, O. Sokolsky, and I. Lee, “Confidence composition for monitors of verification assumptions,” in _2022 ACM/IEEE 13th International Conference on Cyber-Physical Systems (ICCPS)_ , 2022, pp. 1–12.
  * [35] R. G. Newcombe, “Interval estimation for the difference between independent proportions: comparison of eleven methods,” _Statistics in medicine_ , vol. 17, no. 8, pp. 873–890, 1998.

## Appendix A

In this appendix, we provide the proofs of the following theorems from paper
“Causal Repair of Learning-Enabled Cyber-Physical Systems”: Theorem 1 from
Section IV-B, Theorem 2 from Section IV-C, Theorem 3 from Section V-A, and
Theorem 4 from Section V-B.

First, we prove that Algorithm 1 is guaranteed to terminate.

######  Theorem 1 (Termination of Discretization).

Algorithm 1 is guaranteed to terminate.

###### Proof.

Within an input cell $\mathcal{X}_{i}$ of width $\Delta_{x}$, for any two
$x_{1}$, $x_{2}$ inside this cell, we have

$||x_{1}-x_{2}||_{d_{\mathcal{X}}}\leq\sqrt{\text{dim}(\mathcal{X})}\Delta_{x},$
(10)

i.e., any two points within this cell must be smaller than the main diagonal
length of cell $\mathcal{X}_{i}$.

Next, by Assumption 2, the behavior of factual component $f$ is Lipschitz-
continuous, so $\exists$ Lipchitiz constant $c>0$, that

$\forall x_{1},x_{2}\in\mathcal{X},||f(x_{1})-f(x_{2})||_{d_{\mathcal{Y}}}\leq
c||x_{1}-x_{2}||_{d_{\mathcal{X}}}.$ (11)

Therefore, when these inputs are mapped into the output space by $f$, we have

$||f(x_{1})-f(x_{2})||_{d_{\mathcal{Y}}}\leq
c\sqrt{\text{dim}(\mathcal{X})}\Delta_{x},$ (12)

which means there exists a cubic output cell $\mathcal{Y}_{j}$ of width
$\Delta_{y}\geq
c\left(\sqrt{\text{dim}(\mathcal{X})}/\sqrt{\text{dim}(\mathcal{Y})}\right)\Delta_{x}$
such that every input from $\mathcal{X}_{i}$ can be mapped to
$\mathcal{Y}_{j}$.

In this cell partition , the maximal difference between any output
$y\in\mathcal{Y}_{j}$ and its cell center $\text{center}(\mathcal{Y}_{j})$ is
bounded by half of the main diagonal length, i.e.,

$||y-\text{center}(\mathcal{Y}_{j})||_{d_{\mathcal{Y}}}\leq\sqrt{\text{dim}(\mathcal{Y})}\Delta_{y}.$
(13)

Since every output is represented by its cell center, by Definition 2,

$||f-f_{r}||_{d_{f}}\leq\sqrt{\text{dim}(\mathcal{Y})}\Delta_{y}.$ (14)

Recall that Assumption 3 states that the outcome of simulate maintains the
same when two functions has a separation smaller than $\epsilon$. Therefore,
we shrink $\Delta_{x}$ (in lines 7-10) until
$\sqrt{\text{dim}(\mathcal{Y})}\Delta_{y}\leq\epsilon$, so that the
representative function $f_{r}$ leads to the same outcome, i.e. runtime
property violation, as the factual behavior $f$.

Therefore, we are able to find the representative function $f_{r}$ constructed
at line 12. ∎

Next, we prove that the propositional HP model constructed by Algorithm 2
preserves the partial order from Definition 6.

######  Theorem 2 (Encoding Preserves Partial Order).

The encoding $\textsc{encode}:\mathcal{F}_{r}\mapsto\mathcal{V}_{\text{io}}$
specified by evaluating the proposition of every
$u_{ijk}\in\mathcal{U}_{\text{io}}$,as constructed in Algorithm 2, preserves
the partial order as defined in Definition 6.

###### Proof.

We start by proving the forward direction $\Longrightarrow$ in Definition 6.
Suppose we have two arbitrary behavior choices choices
$f_{1},f_{2}\in\mathcal{F}_{r}$. The factual behavior is represented by
$f_{r}$ and we have $f_{1}\preccurlyeq_{f_{r}}f_{2}$.

By Definition 4, this partial order $\preccurlyeq_{f_{r}}$ means that $f_{1}$
has an output in between that of $f_{r}$ and $f_{2}$ in every output
dimension, i.e.,

$\begin{split}&\forall x\in\mathcal{X},\forall
j\in\\{1,\dots,\text{dim}(\mathcal{Y})\\},\\\ &(f_{r}(x)[j]\leq
f_{1}(x)[j]\leq f_{2}(x)[j])\\\ &\lor(f_{r}(x)[j]\geq f_{1}(x)[j]\geq
f_{2}(x)[j])\end{split}$ (15)

We pick an arbitrary input $x\in$ input cell $\mathcal{X}_{i}$, and an
arbitrary output dimension $j$. By Equation 15, we know that the output bins
that $f_{r}(x)$, $f_{1}(x)$ and $f_{2}(x)$ in dimension $j$ must be in either
non-decreasing, or non-increasing order. Without loss of generality, we
consider non-decreasing order here. That is,

$\begin{split}&f_{r}(x)\in bin_{\mathcal{Y}}(j,k_{r})\\\ &f_{1}(x)\in
bin_{\mathcal{Y}}(j,k_{1})\\\ &f_{2}(x)\in bin_{\mathcal{Y}}(j,k_{2})\\\
&k_{r}\leq k_{1}\leq k_{3}\end{split}$ (16)

Therefore, we have $lo(bin_{\mathcal{Y}}(j,k_{r}))\leq
lo(bin_{\mathcal{Y}}(j,k_{1}))\leq lo(bin_{\mathcal{Y}}(j,k_{2}))$. Now we
consider the encoded node values on the subset of nodes
$\mathcal{U}_{\text{io}}(i,j)\subseteq\mathcal{U}_{\text{io}}$, which are the
nodes that encodes the behaviors of input cell $\mathcal{X}_{i}$ on dimension
$j$, i.e., $\mathcal{U}_{\text{io}}(i,j)=\\{u_{ijk}\mid k=1,\dots,n_{j}\\}$.
The value assignments on $\mathcal{U}_{\text{io}}(i,j)$ given by the three
behavior choices are denoted as $\mathbf{v}_{r}(i,j)$, $\mathbf{v}_{1}(i,j)$
and $\mathbf{v}_{2}(i,j)\in\mathcal{V}_{\text{io}}(i,j)$, respectively.

Based on the order $lo(bin_{\mathcal{Y}}(j,k_{r}))\leq
lo(bin_{\mathcal{Y}}(j,k_{1}))\leq lo(bin_{\mathcal{Y}}(j,k_{2}))$, evaluating
the node propositions, i.e., the encoding function $\textsc{encode}()$,
results in

$\begin{split}&\mathbf{v}_{r}(i,j)=\underbrace{111...1}_{k_{r}}\underbrace{000...0}_{n_{j}-k_{r}}\\\
&\mathbf{v}_{1}(i,j)=\underbrace{111...11}_{k_{1}}\underbrace{00...0}_{n_{j}-k_{1}}\\\
&\mathbf{v}_{2}(i,j)=\underbrace{111...111}_{k_{2}}\underbrace{0...0}_{n_{j}-k_{2}}\end{split}$
(17)

From the evaluated assignments, we can see the two susbets of differing nodes
between $\mathbf{v}_{r}$ and $\mathbf{v}_{1}$ and between $\mathbf{v}_{r}$ and
$\mathbf{v}_{2}$ are

$\begin{split}\text{diff}_{\mathcal{U}_{\text{io}}}(i,j)(\mathbf{v}_{r},\mathbf{v}_{1})=\\{u_{ijk}\mid
k_{r}+1\leq k\leq k_{1}\\}\\\
\text{diff}_{\mathcal{U}_{\text{io}}}(i,j)(\mathbf{v}_{r},\mathbf{v}_{2})=\\{u_{ijk}\mid
k_{r}+1\leq k\leq k_{2}\\}\end{split}$ (18)

And we can see

$\text{diff}_{\mathcal{U}_{\text{io}}}(i,j)(\mathbf{v}_{r},\mathbf{v}_{1})\subseteq\text{diff}_{\mathcal{U}_{\text{io}}}(i,j)(\mathbf{v}_{r},\mathbf{v}_{2})$
(19)

The same set containment holds when we use $\geq$, i.e. non-increasing order.

Notice that the three integers $k_{r}$, $k_{1}$ and $k_{2}$ denote the number
of nodes flipped to 1 in $\mathcal{U}_{\text{io}}(i,j)$ by these 3 behavior
choices. Therefore, we can as well denote them as $k_{r}(i,j)$, $k_{1}(i,j)$
and $k_{2}(i,j)$. Since Equation 15 applies for all input cells
$\mathcal{X}_{i}$ and output dimensions $j$, we have $k_{r}(i,j)\leq
k_{1}(i,j)\leq k_{2}(i,j)$ (or $\geq$) for each pair of $i$ and $j$.
Therefore, on all $i,j$,
$\text{diff}_{\mathcal{U}_{\text{io}}}(i,j)(\mathbf{v}_{r},\mathbf{v}_{1})\subseteq\text{diff}_{\mathcal{U}_{\text{io}}}(i,j)(\mathbf{v}_{r},\mathbf{v}_{2})$
holds. Consequently,

$\begin{split}&\text{diff}_{\mathcal{U}_{\text{io}}}(\mathbf{v}_{r},\mathbf{v}_{1})=\bigcup_{i,j}\text{diff}_{\mathcal{U}_{\text{io}}}(i,j)(\mathbf{v}_{r},\mathbf{v}_{1})\\\
\subseteq&\bigcup_{i,j}\text{diff}_{\mathcal{U}_{\text{io}}}(i,j)(\mathbf{v}_{r},\mathbf{v}_{2})=\text{diff}_{\mathcal{U}_{\text{io}}}(\mathbf{v}_{r},\mathbf{v}_{2})\end{split}$
(20)

By Definition 5, we have
$\mathbf{v}_{1}\preccurlyeq_{\mathbf{v}_{r}}\mathbf{v}_{2}$, and the forward
direction $\Longrightarrow$ is proved.

To prove the backward direction $\Longleftarrow$, we can use a symmetrical
approach as above. That is, we still split $\mathcal{U}_{\text{io}}$ into
multiple $\mathcal{U}_{\text{io}}(i,j)$. Since we assume
$\text{diff}_{\mathcal{U}_{\text{io}}}(\mathbf{v}_{r},\mathbf{v}_{1})\subseteq\text{diff}_{\mathcal{U}_{\text{io}}}(\mathbf{v}_{r},\mathbf{v}_{2})$,
we must have $\forall
i,j,\text{diff}_{\mathcal{U}_{\text{io}}}(i,j)(\mathbf{v}_{r},\mathbf{v}_{1})\subseteq\text{diff}_{\mathcal{U}_{\text{io}}}(i,j)(\mathbf{v}_{r},\mathbf{v}_{2})$.
Consequently, we can show that $f_{1}(x)[j]\leq f_{2}(x)[j]$ (or $\geq$) on
arbitrary input space $x\in\mathcal{X}_{i}$ and output dimension $j$ based on
the construction in Algorithm 2, and therefore
$f_{1}\preccurlyeq_{f_{r}}f_{2}$. This encoding indeed preserves partial order
in both directions.

∎

Next, we prove the validity of the probabilistic statement reported by
Algorithm 3 when it cannot find a suitable counter example.

######  Theorem 3 (Probabilistic Guarantee on Search Failure).

Given an HP model constructed by Algorithm 2, a probability threshold
$p\in[0,1]$ and a confidence $1-\alpha\in[0,1]$, the final statement at line 9
of Algorithm 3 holds.

###### Proof.

Let $\mathcal{D}_{1}$ denote the uniform distribution on value assignments on
$\mathcal{U}_{\text{io}}$, i.e., all propositional assignments on the binary
nodes. Consider random variable $\mathbf{v}\sim\mathcal{D}_{1}$. Let
$s=\textsc{simulate}(\textsc{decode}(\mathbf{v}))$ denote whether
$\mathcal{S}$ is safe or not when component $\mathcal{C}$ is replaced by a
counter factual $f^{\prime}=\textsc{decode}(\mathbf{v})$. Note that $s$ is a
Bernoulli random variable, as it has value $1$ when $\mathcal{S}$ is safe and
$0$ otherwise and that $Pr_{\mathcal{D}_{1}}[s=1]$ is the number of safe
assignments divided by the total number of feasible assignments $n^{m}$.

For a Bernoulli distribution, the number of successes, denoted $N_{s}$, from
$N$ sequential samples follows a Binomial distribution. If $N_{s}=0$ we can
apply the Wilson score interval [35] to get an upper bound on the success
probability of the Bernoulli random variable:

$\begin{split}z&=Q(1-\frac{\alpha}{2})\\\
\Pr[\text{success}]&\in\frac{N_{s}+\frac{z^{2}}{2}}{N+z^{2}}\pm\frac{z}{N+z^{2}}\sqrt{\frac{N_{s}(N-N_{s})}{N}+\frac{z^{2}}{4}}\\\
&=\frac{z^{2}}{2N+2z^{2}}\pm\frac{z^{2}}{2N+2z^{2}}\\\
&=[0,\frac{z^{2}}{N+z^{2}}]\end{split}$ (21)

Equation (21) means that, if we do not see any success in $N$ Bernoulli
trials, then with $1-\alpha$ confidence we estimate that the probability of a
success in one trial is less than $z^{2}/(N+z^{2})$, with $z$ being the
$(1-\alpha/2)$-quantile of standard normal distribution.

Therefore, given an upper-bound probability estimation $p$, we can compute $N$
in terms of $p$ and $\alpha$.

$\begin{split}&p=\frac{z^{2}}{(N+z^{2})}\\\
\implies&N=(\frac{1}{p}-1)z^{2}=(\frac{1}{p}-1)Q(1-\frac{\alpha}{2})^{2}\end{split}$
(22)

Therefore, if we sample at least $(1/p-1)Q(1-\alpha/2)^{2}$ value assignments
without any success, we can make the claim in line 9, supported by Wilson
score interval. Notice that this minimal number of sampling is met at line 1.
Consequently, the statement at line 9 holds. ∎

In the end, we prove the output of Algorithm 4 is an actual cause given HP
model $\mathcal{M}$.

######  Theorem 4 (Output is Actual Cause).

Let HP model $\mathcal{M}$ constructed in Algorithm 2 be given with factual
node value assignment $\mathbf{v}_{\text{io}}$. Let
$\mathbf{v}_{\text{io}}^{\prime}$ be a counterfactual node value assignment
from Algorithm 3. The node values on $\mathcal{U}_{\text{io}}$ where
assignments $\mathbf{v}_{\text{io}}$ and $\mathbf{v}_{\text{io}}^{*}$ disagree
in Algorithm 4 are an _actual cause_ of $\mathcal{S}(f)\not\models\varphi$ as
per $\mathcal{M}$ constructed in Algorithm 2.

###### Proof.

First, the encoding of the factual function, $\mathbf{v}_{\text{io}}$ gives
property violation $\neg\varphi$ based on Algorithm 1 and 2. Therefore, AC1
holds.

Second, if we partition $\mathcal{U}_{\text{endo}}$ into
$\mathcal{U}_{1}=\operatorname{diff}_{\mathcal{U}_{\text{io}}}^{*}\cup\\{u_{\varphi}\\}$,
$\mathcal{U}_{2}=\emptyset$ and
$\mathcal{U}_{3}=\mathcal{U}_{\text{endo}}\setminus\mathcal{U}_{1}$, we can
check this partition on AC2. The subset $\mathcal{U}_{3}$ does not have
counterfactual values in these assignments so there is no difference between
the counterfactuals in AC2(a) and AC2(b). In fact, the existence of this
partition means AC2 holds.

Finally, there does not exist a subset
$\text{diff}_{\mathcal{U}_{\text{io}}}^{**}\subset\text{diff}_{\mathcal{U}_{\text{io}}}^{*}$
such that by using counterfactual values on this smaller subset only gives
satisfaction of $\varphi$, as then Algorithm 3 would have output the node
value assignment corresponding to that counterfactual instead of outputting
$\mathbf{v}_{\text{io}}^{\prime}$.

Since all the three conditions hold, the factual value assignments on the
disagreeing nodes, i.e., $\text{diff}_{\mathcal{U}_{\text{io}}}^{*}$, is an
actual cause of the property violation. ∎