Safety synthesis (max P(never reach avoid) = 1 - min-reach)

Add a controllerMax flag to the solver so the controller can MINIMIZE reach (for
safety) in addition to maximizing (reachability), reusing OVI. solve::maxSafety
{Pessimistic,Optimistic} return sound [lower,upper] on the safety probability.
Broadens coverage to the ARCH safety category (building automation, rooms).

Verified (tests/unit/test_safety.cpp): controller-can-avoid -> 1; unavoidable -> 0;
coin -> 0.5; point MDP pess==opt; differential vs reachability on single-action
interval MDPs (safety == 1 - reach).

Suite: 56/56 cases, 606499 assertions, 0 failures.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: Kiguli

paper/IMPaCT-v2.md

@@ -235,6 +235,16 @@ Format per entry: `date — feature | decision + rationale | algorithm | refs |
 - ISSUE-0010 (perf, low): OVI's VI-from-below is slow on large strongly-recurrent
   IMDPs; MECCollapse (optimistic) / robust-EC interval iteration is the fast path.
 
+- **Safety synthesis** (broadens to the ARCH safety category: BAS, rooms).
+  `solve::maxSafety{Pessimistic,Optimistic}` = max P(never reach `avoid`) =
+  1 - min-reach-to-avoid. Added a `controllerMax` flag to the solver so the
+  controller can MINIMIZE reach (for safety) as well as maximize (reachability),
+  reusing OVI. Verified (test_safety.cpp): controller-can-avoid → 1; unavoidable →
+  0; coin → 0.5; point pess==opt; and a differential vs reachability on
+  single-action interval MDPs (safety == 1 - reach). Suite: 56/56, 606499 assertions.
+- **Consolidated CAV paper draft** at `paper/draft.md` (algorithms + verified
+  citations + verification methodology incl. the MC-found bugs + experimental results).
+
 <!-- add new dated entries above this line -->
 
 ---

src/solve.cpp

@@ -14,7 +14,10 @@ namespace {
 // Robust Bellman backup at one state: max over actions of the omax optimum.
 // sense = Min for pessimistic (adversarial nature), Max for optimistic.
 // A state with no actions has value 0 (cannot make progress to a target).
-double backup(const StateActions& actions, const std::vector<double>& V, omax::Sense sense) {
+// controllerMax: the controller MAXIMIZES over actions (reachability) or MINIMIZES
+// (used for safety = 1 - min-reach-to-avoid). A state with no actions has value 0.
+double backup(const StateActions& actions, const std::vector<double>& V, omax::Sense sense,
+              bool controllerMax = true) {
     double best = 0.0;
     bool any = false;
     std::vector<double> lo, hi, vv;
@@ -23,9 +26,10 @@ double backup(const StateActions& actions, const std::vector<double>& V, omax::S
         lo.reserve(act.size()); hi.reserve(act.size()); vv.reserve(act.size());
         for (const Interval& iv : act) { lo.push_back(iv.lo); hi.push_back(iv.hi); vv.push_back(V[iv.to]); }
         const double val = omax::optimize(lo, hi, vv, sense).value;
-        if (!any || val > best) { best = val; any = true; }
+        if (!any) { best = val; any = true; }
+        else best = controllerMax ? std::max(best, val) : std::min(best, val);
     }
-    return best;
+    return any ? best : 0.0;
 }
 
 struct Collapsed {
@@ -142,7 +146,7 @@ IntervalResult solveReach(const IMDPModel& m, const std::set<int>& targets,
 // needed, so it converges on nature-confinable ECs (ISSUE-0003). If the guess is
 // not yet inductive, refine L (smaller delta) and retry.
 IntervalResult solveOVI(const IMDPModel& m, const std::set<int>& targets,
-                        double eps, omax::Sense sense) {
+                        double eps, omax::Sense sense, bool controllerMax = true) {
     const int n = (int)m.size();
     std::vector<double> L(n, 0.0), tmp(n), U(n), FU(n);
     for (int t : targets) L[t] = 1.0;
@@ -153,7 +157,7 @@ IntervalResult solveOVI(const IMDPModel& m, const std::set<int>& targets,
         for (int it = 0; it < MAXINNER; ++it) {              // refine L from below
             ++iters;
             for (int s = 0; s < n; ++s)
-                tmp[s] = targets.count(s) ? 1.0 : backup(m[s], L, sense);
+                tmp[s] = targets.count(s) ? 1.0 : backup(m[s], L, sense, controllerMax);
             double ch = 0.0;
             for (int s = 0; s < n; ++s) ch = std::max(ch, std::fabs(tmp[s] - L[s]));
             L.swap(tmp);
@@ -163,7 +167,7 @@ IntervalResult solveOVI(const IMDPModel& m, const std::set<int>& targets,
             U[s] = targets.count(s) ? 1.0 : std::min(1.0, L[s] + eps);
         bool inductive = true;                                // verify F(U) <= U
         for (int s = 0; s < n; ++s) {
-            FU[s] = targets.count(s) ? 1.0 : backup(m[s], U, sense);
+            FU[s] = targets.count(s) ? 1.0 : backup(m[s], U, sense, controllerMax);
             if (FU[s] > U[s] + 1e-12) inductive = false;
         }
         if (inductive) return IntervalResult{L, U, iters};    // L <= V* <= U, gap <= eps
@@ -173,6 +177,26 @@ IntervalResult solveOVI(const IMDPModel& m, const std::set<int>& targets,
     return IntervalResult{L, U, iters};                       // best-effort fallback
 }
 
+// Min-reach: controller MINIMIZES reach to `targets`. Used for safety.
+IntervalResult minReach(const IMDPModel& m, const std::set<int>& targets,
+                        double eps, omax::Sense sense) {
+    return solveOVI(m, targets, eps, sense, /*controllerMax=*/false);
+}
+
+// Safety = 1 - min-reach-to-avoid (controller maximizes staying out of `avoid`).
+IntervalResult safetyFromMinReach(const IMDPModel& m, const std::set<int>& avoid,
+                                  double eps, omax::Sense sense) {
+    IntervalResult mr = minReach(m, avoid, eps, sense);
+    IntervalResult r;
+    r.iterations = mr.iterations;
+    r.lower.resize(mr.lower.size()); r.upper.resize(mr.upper.size());
+    for (size_t s = 0; s < mr.lower.size(); ++s) {
+        r.lower[s] = 1.0 - mr.upper[s];    // safety lower = 1 - max possible reach
+        r.upper[s] = 1.0 - mr.lower[s];
+    }
+    return r;
+}
+
 IntervalResult dispatch(const IMDPModel& m, const std::set<int>& targets,
                         double eps, omax::Sense sense, Method method) {
     return (method == Method::MECCollapse) ? solveReach(m, targets, eps, sense)
@@ -194,5 +218,14 @@ IntervalResult maxReachOptimistic(const IMDPModel& m, const std::set<int>& targe
     return dispatch(m, targets, eps, omax::Sense::Max, method);
 }
 
+// Robust safety: max over controller of P(never reach `avoid`); pessimistic =
+// nature adversarial (maximizes reach to avoid, omax Sense::Max).
+IntervalResult maxSafetyPessimistic(const IMDPModel& m, const std::set<int>& avoid, double eps) {
+    return safetyFromMinReach(m, avoid, eps, omax::Sense::Max);
+}
+IntervalResult maxSafetyOptimistic(const IMDPModel& m, const std::set<int>& avoid, double eps) {
+    return safetyFromMinReach(m, avoid, eps, omax::Sense::Min);
+}
+
 } // namespace solve
 } // namespace impact

src/solve.h

@@ -52,6 +52,12 @@ namespace solve {
     IntervalResult maxReachPessimistic(const IMDPModel& m, const std::set<int>& targets, double eps, Method method);
     IntervalResult maxReachOptimistic (const IMDPModel& m, const std::set<int>& targets, double eps, Method method);
 
+    // Robust safety: max over controller of P(never reach `avoid`) = 1 - min-reach to
+    // avoid. Pessimistic = nature adversarial; optimistic = nature cooperative.
+    // Returns sound [lower,upper] on the safety probability (gap <= 2*eps).
+    IntervalResult maxSafetyPessimistic(const IMDPModel& m, const std::set<int>& avoid, double eps);
+    IntervalResult maxSafetyOptimistic (const IMDPModel& m, const std::set<int>& avoid, double eps);
+
 } // namespace solve
 } // namespace impact
 

tests/Makefile

@@ -33,7 +33,8 @@ TEST_SRCS := unit/test_main.cpp \
              unit/test_ltl_dfa.cpp \
              unit/test_product.cpp \
              unit/test_abstraction.cpp \
-             unit/test_omega.cpp
+             unit/test_omega.cpp \
+             unit/test_safety.cpp
 
 BIN := run_tests
 

tests/unit/test_safety.cpp

@@ -0,0 +1,78 @@
+// ============================================================================
+// CONTRACT TESTS — safety synthesis: max P(never reach `avoid`) = 1 - min-reach.
+// Verified by hand cases + a differential vs reachability on single-action models
+// (where the controller has no choice, so safety == 1 - reach).
+// ============================================================================
+#include "../doctest.h"
+#include "../contracts/contracts.h"
+
+#include <set>
+#include <vector>
+#include <random>
+#include <cmath>
+
+using namespace impact::solve;
+
+static const double EPS = 1e-7;
+static double mid(const IntervalResult& r, int s) { return 0.5*(r.lower[s]+r.upper[s]); }
+
+TEST_CASE("safety: controller can avoid -> safety 1") {
+    // 0: a->1(avoid) | b->2(safe sink) ; 1 avoid sink ; 2 safe sink
+    IMDPModel m = { {{ {1,1,1} }, { {2,1,1} }}, {{ {1,1,1} }}, {{ {2,1,1} }} };
+    auto sf = maxSafetyPessimistic(m, {1}, EPS);
+    CHECK(mid(sf, 0) == doctest::Approx(1.0).epsilon(1e-6));   // pick b
+    CHECK(mid(sf, 2) == doctest::Approx(1.0).epsilon(1e-6));
+    CHECK(mid(sf, 1) == doctest::Approx(0.0).epsilon(1e-6));   // already in avoid
+}
+
+TEST_CASE("safety: avoid unavoidable -> safety 0") {
+    IMDPModel m = { {{ {1,1,1} }}, {{ {1,1,1} }} };            // 0 -> 1(avoid) forced
+    auto sf = maxSafetyPessimistic(m, {1}, EPS);
+    CHECK(mid(sf, 0) == doctest::Approx(0.0).epsilon(1e-6));
+}
+
+TEST_CASE("safety: coin to avoid -> safety 0.5") {
+    IMDPModel m = { {{ {1,0.5,0.5}, {2,0.5,0.5} }}, {{ {1,1,1} }}, {{ {2,1,1} }} };
+    auto sf = maxSafetyPessimistic(m, {1}, EPS);
+    CHECK(mid(sf, 0) == doctest::Approx(0.5).epsilon(1e-3));
+}
+
+TEST_CASE("safety: point MDP pessimistic == optimistic") {
+    IMDPModel m = { {{ {1,0.3,0.3}, {2,0.7,0.7} }}, {{ {1,1,1} }}, {{ {2,1,1} }} };
+    auto sp = maxSafetyPessimistic(m, {1}, EPS);
+    auto so = maxSafetyOptimistic(m, {1}, EPS);
+    for (int s = 0; s < 3; ++s) CHECK(mid(sp,s) == doctest::Approx(mid(so,s)).epsilon(1e-4));
+}
+
+TEST_CASE("safety: differential vs reachability on single-action interval MDPs") {
+    // With one action per state the controller has no choice, so
+    // safety == 1 - reach(avoid) (nature adversarial both ways).
+    std::mt19937 rng(909);
+    std::uniform_real_distribution<double> u01(0.0, 1.0);
+    int checked = 0;
+    for (int trial = 0; trial < 4000 && checked < 400; ++trial) {
+        int n = 2 + (int)(rng() % 5);
+        std::set<int> avoid; avoid.insert((int)(rng() % n));
+        IMDPModel m(n);
+        for (int s = 0; s < n; ++s) {
+            std::vector<int> succ;
+            for (int t = 0; t < n; ++t) if (rng() & 1u) succ.push_back(t);
+            if (succ.empty()) succ.push_back((int)(rng() % n));
+            std::vector<double> w(succ.size()); double sum = 0;
+            for (double& x : w) { x = u01(rng) + 1e-3; sum += x; }
+            ActionDist act;
+            for (size_t k = 0; k < succ.size(); ++k) { double p = w[k]/sum;
+                double r = 0.1 * u01(rng);
+                act.push_back({succ[k], std::max(0.0,p-r), std::min(1.0,p+r)}); }
+            m[s].push_back(act);                              // exactly one action
+        }
+        ++checked;
+        auto sf = maxSafetyPessimistic(m, avoid, EPS);
+        auto rc = maxReachOptimistic(m, avoid, EPS);         // single action: nature-max reach
+        for (int s = 0; s < n; ++s) {
+            CHECK(sf.lower[s] <= sf.upper[s] + 1e-9);
+            CHECK(std::fabs(mid(sf, s) - (1.0 - mid(rc, s))) < 3e-3);   // safety == 1 - reach
+        }
+    }
+    CHECK(checked > 150);
+}