Theory HS_VC_Examples

(*  Title:       Examples of hybrid systems verifications
    Author:      Jonathan Julián Huerta y Munive, 2020
    Maintainer:  Jonathan Julián Huerta y Munive <jonjulian23@gmail.com>
*)

subsection ‹ Examples ›

text ‹ We prove partial correctness specifications of some hybrid systems with our
verification components.›

theory HS_VC_Examples
  imports HS_VC_Spartan

begin


subsubsection ‹Pendulum›

text ‹ The ODEs @{text "x' t = y t"} and {text "y' t = - x t"} describe the circular motion of
a mass attached to a string looked from above. We use @{text "s$1"} to represent the x-coordinate
and @{text "s$2"} for the y-coordinate. We prove that this motion remains circular. ›

abbreviation fpend :: "real^2  real^2" ("f")
  where "f s  (χ i. if i = 1 then s$2 else -s$1)"

abbreviation pend_flow :: "real  real^2  real^2" ("φ")
  where "φ t s  (χ i. if i = 1 then s$1 * cos t + s$2 * sin t else - s$1 * sin t + s$2 * cos t)"

― ‹Verified with annotated dynamics. ›

lemma pendulum_dyn: "(λs. r2 = (s$1)2 + (s$2)2)  |EVOL φ G T] (λs. r2 = (s$1)2 + (s$2)2)"
  by force

― ‹Verified with differential invariants. ›

lemma pendulum_inv: "(λs. r2 = (s$1)2 + (s$2)2)  |x´= f & G] (λs. r2 = (s$1)2 + (s$2)2)"
  by (auto intro!: diff_invariant_rules poly_derivatives)

― ‹Verified with the flow. ›

lemma local_flow_pend: "local_flow f UNIV UNIV φ"
  apply(unfold_locales, simp_all add: local_lipschitz_def lipschitz_on_def vec_eq_iff, clarsimp)
    apply(rule_tac x="1" in exI, clarsimp, rule_tac x=1 in exI)
    apply(simp add: dist_norm norm_vec_def L2_set_def power2_commute UNIV_2)
  by (auto simp: forall_2 intro!: poly_derivatives)

lemma pendulum_flow: "(λs. r2 = (s$1)2 + (s$2)2)  |x´=f & G] (λs. r2 = (s$1)2 + (s$2)2)"
  by (force simp: local_flow.fbox_g_ode_subset[OF local_flow_pend])

no_notation fpend ("f")
        and pend_flow ("φ")


subsubsection ‹ Bouncing Ball ›

text ‹ A ball is dropped from rest at an initial height @{text "h"}. The motion is described with
the free-fall equations @{text "x' t = v t"} and @{text "v' t = g"} where @{text "g"} is the
constant acceleration due to gravity. The bounce is modelled with a variable assignment that
flips the velocity. That is, we model it as a completely elastic collision with the ground. We use 
@{text "s$1"} to represent the ball's height and @{text "s$2"} for its velocity. We prove that the 
ball remains above ground and below its initial resting position. ›

abbreviation fball :: "real  real^2  real^2" ("f")
  where "f g s  (χ i. if i = 1 then s$2 else g)"

abbreviation ball_flow :: "real  real  real^2  real^2" ("φ")
  where "φ g t s  (χ i. if i = 1 then g * t ^ 2/2 + s$2 * t + s$1 else g * t + s$2)"

― ‹Verified with differential invariants. ›

named_theorems bb_real_arith "real arithmetic properties for the bouncing ball."

lemma inv_imp_pos_le[bb_real_arith]:
  assumes "0 > g" and inv: "2 * g * x - 2 * g * h = v * v"
  shows "(x::real)  h"
proof-
  have "v * v = 2 * g * x - 2 * g * h  0 > g"
    using inv and 0 > g by auto
  hence obs:"v * v = 2 * g * (x - h)  0 > g  v * v  0"
    using left_diff_distrib mult.commute by (metis zero_le_square)
  hence "(v * v)/(2 * g) = (x - h)"
    by auto
  also from obs have "(v * v)/(2 * g)  0"
    using divide_nonneg_neg by fastforce
  ultimately have "h - x  0"
    by linarith
  thus ?thesis by auto
qed

lemma "diff_invariant (λs. 2 * g * s$1 - 2 * g * h - s$2 * s$2 = 0) (λt. f g) (λs. UNIV) S t0 G"
  by (auto intro!: poly_derivatives diff_invariant_rules)

lemma bouncing_ball_inv: "g < 0  h  0 
  (λs. s$1 = h  s$2 = 0) 
  |LOOP (
    (x´=(f g) & (λ s. s$1  0) DINV (λs. 2 * g * s$1 - 2 * g * h - s$2 * s$2 = 0)) ;
    (IF (λ s. s$1 = 0) THEN (2 ::= (λs. - s$2)) ELSE skip))
  INV (λs. 0  s$1  2 * g * s$1 - 2 * g * h - s$2 * s$2 = 0)]
  (λs. 0  s$1  s$1  h)"
  apply(rule fbox_loopI, simp_all, force, force simp: bb_real_arith)
  by (rule fbox_g_odei) (auto intro!: poly_derivatives diff_invariant_rules)

― ‹Verified with annotated dynamics. ›

lemma inv_conserv_at_ground[bb_real_arith]:
  assumes invar: "2 * g * x = 2 * g * h + v * v"
    and pos: "g * τ2 / 2 + v * τ + (x::real) = 0"
  shows "2 * g * h + (g * τ + v) * (g * τ + v) = 0"
proof-
  from pos have "g * τ2  + 2 * v * τ + 2 * x = 0" by auto
  then have "g2 * τ2  + 2 * g * v * τ + 2 * g * x = 0"
    by (metis (mono_tags) Groups.mult_ac(1,3) mult_zero_right
        monoid_mult_class.power2_eq_square semiring_class.distrib_left)
  hence "g2 * τ2 + 2 * g * v * τ + v2 + 2 * g * h = 0"
    using invar by (simp add: monoid_mult_class.power2_eq_square)
  hence obs: "(g * τ + v)2 + 2 * g * h = 0"
    apply(subst power2_sum) by (metis (no_types) Groups.add_ac(2, 3)
        Groups.mult_ac(2, 3) monoid_mult_class.power2_eq_square nat_distrib(2))
  thus "2 * g * h + (g * τ + v) * (g * τ + v) = 0"
    by (simp add: add.commute distrib_right power2_eq_square)
qed

lemma inv_conserv_at_air[bb_real_arith]:
  assumes invar: "2 * g * x = 2 * g * h + v * v"
  shows "2 * g * (g * τ2 / 2 + v * τ + (x::real)) =
  2 * g * h + (g * τ + v) * (g * τ + v)" (is "?lhs = ?rhs")
proof-
  have "?lhs = g2 * τ2 + 2 * g * v * τ + 2 * g * x"
    by(auto simp: algebra_simps semiring_normalization_rules(29))
  also have "... = g2 * τ2 + 2 * g * v * τ + 2 * g * h + v * v" (is "... = ?middle")
    by(subst invar, simp)
  finally have "?lhs = ?middle".
  moreover
  {have "?rhs = g * g * (τ * τ) + 2 * g * v * τ + 2 * g * h + v * v"
    by (simp add: Groups.mult_ac(2,3) semiring_class.distrib_left)
  also have "... = ?middle"
    by (simp add: semiring_normalization_rules(29))
  finally have "?rhs = ?middle".}
  ultimately show ?thesis by auto
qed

lemma bouncing_ball_dyn: "g < 0  h  0 
  (λs. s$1 = h  s$2 = 0) 
  |LOOP (
    (EVOL (φ g) (λ s. s$1  0) T) ;
    (IF (λ s. s$1 = 0) THEN (2 ::= (λs. - s$2)) ELSE skip))
  INV (λs. 0  s$1 2 * g * s$1 = 2 * g * h + s$2 * s$2)]
  (λs. 0  s$1  s$1  h)"
  by (rule fbox_loopI) (auto simp: bb_real_arith)

― ‹Verified with the flow. ›

lemma local_flow_ball: "local_flow (f g) UNIV UNIV (φ g)"
  apply(unfold_locales, simp_all add: local_lipschitz_def lipschitz_on_def vec_eq_iff, clarsimp)
    apply(rule_tac x="1/2" in exI, clarsimp, rule_tac x=1 in exI)
    apply(simp add: dist_norm norm_vec_def L2_set_def UNIV_2)
  by (auto simp: forall_2 intro!: poly_derivatives)

lemma bouncing_ball_flow: "g < 0  h  0 
  (λs. s$1 = h  s$2 = 0) 
  |LOOP (
    (x´=(λt. f g) & (λ s. s$1  0) on (λs. UNIV) UNIV @ 0) ;
    (IF (λ s. s$1 = 0) THEN (2 ::= (λs. - s$2)) ELSE skip))
  INV (λs. 0  s$1 2 * g * s$1 = 2 * g * h + s$2 * s$2)]
  (λs. 0  s$1  s$1  h)"
  apply(rule fbox_loopI, simp_all add: local_flow.fbox_g_ode_subset[OF local_flow_ball])
  by (auto simp: bb_real_arith)

no_notation fball ("f")
        and ball_flow ("φ")


subsubsection ‹ Thermostat ›

text ‹ A thermostat has a chronometer, a thermometer and a switch to turn on and off a heater.
At most every @{text "t"} minutes, it sets its chronometer to @{term "0::real"}, it registers
the room temperature, and it turns the heater on (or off) based on this reading. The temperature
follows the ODE @{text "T' = - a * (T - U)"} where @{text "U"} is @{text "L ≥ 0"} when the heater
is on, and @{text "0"} when it is off. We use @{term "1::4"} to denote the room's temperature,
@{term "2::4"} is time as measured by the thermostat's chronometer, @{term "3::4"} is the
temperature detected by the thermometer, and @{term "4::4"} states whether the heater is on
(@{text "s$4 = 1"}) or off (@{text "s$4 = 0"}). We prove that the thermostat keeps the room's
temperature between @{text "Tmin"} and @{text "Tmax"}. ›

abbreviation temp_vec_field :: "real  real  real^4  real^4" ("f")
  where "f a L s  (χ i. if i = 2 then 1 else (if i = 1 then - a * (s$1 - L) else 0))"

abbreviation temp_flow :: "real  real  real  real^4  real^4" ("φ")
  where "φ a L t s  (χ i. if i = 1 then - exp(-a * t) * (L - s$1) + L else
  (if i = 2 then t + s$2 else s$i))"

― ‹Verified with the flow. ›

lemma norm_diff_temp_dyn: "0 < a  f a L s1 - f a L s2 = ¦a¦ * ¦s1$1 - s2$1¦"
proof(simp add: norm_vec_def L2_set_def, unfold UNIV_4, simp)
  assume a1: "0 < a"
  have f2: "r ra. ¦(r::real) + - ra¦ = ¦ra + - r¦"
    by (metis abs_minus_commute minus_real_def)
  have "r ra rb. (r::real) * ra + - (r * rb) = r * (ra + - rb)"
    by (metis minus_real_def right_diff_distrib)
  hence "¦a * (s1$1 + - L) + - (a * (s2$1 + - L))¦ = a * ¦s1$1 + - s2$1¦"
    using a1 by (simp add: abs_mult)
  thus "¦a * (s2$1 - L) - a * (s1$1 - L)¦ = a * ¦s1$1 - s2$1¦"
    using f2 minus_real_def by presburger
qed

lemma local_lipschitz_temp_dyn:
  assumes "0 < (a::real)"
  shows "local_lipschitz UNIV UNIV (λt::real. f a L)"
  apply(unfold local_lipschitz_def lipschitz_on_def dist_norm)
  apply(clarsimp, rule_tac x=1 in exI, clarsimp, rule_tac x=a in exI)
  using assms
  apply(simp add: norm_diff_temp_dyn)
  apply(simp add: norm_vec_def L2_set_def, unfold UNIV_4, clarsimp)
  unfolding real_sqrt_abs[symmetric] by (rule real_le_lsqrt) auto

lemma local_flow_temp: "a > 0  local_flow (f a L) UNIV UNIV (φ a L)"
  by (unfold_locales, auto intro!: poly_derivatives local_lipschitz_temp_dyn simp: forall_4 vec_eq_iff)

lemma temp_dyn_down_real_arith:
  assumes "a > 0" and Thyps: "0 < Tmin" "Tmin  T" "T  Tmax"
    and thyps: "0  (t::real)" "τ{0..t}. τ  - (ln (Tmin / T) / a) "
  shows "Tmin  exp (-a * t) * T" and "exp (-a * t) * T  Tmax"
proof-
  have "0  t  t  - (ln (Tmin / T) / a)"
    using thyps by auto
  hence "ln (Tmin / T)  - a * t  - a * t  0"
    using assms(1) divide_le_cancel by fastforce
  also have "Tmin / T > 0"
    using Thyps by auto
  ultimately have obs: "Tmin / T  exp (-a * t)" "exp (-a * t)  1"
    using exp_ln exp_le_one_iff by (metis exp_less_cancel_iff not_less, simp)
  thus "Tmin  exp (-a * t) * T"
    using Thyps by (simp add: pos_divide_le_eq)
  show "exp (-a * t) * T  Tmax"
    using Thyps mult_left_le_one_le[OF _ exp_ge_zero obs(2), of T]
      less_eq_real_def order_trans_rules(23) by blast
qed

lemma temp_dyn_up_real_arith:
  assumes "a > 0" and Thyps: "Tmin  T" "T  Tmax" "Tmax < (L::real)"
    and thyps: "0  t" "τ{0..t}. τ  - (ln ((L - Tmax) / (L - T)) / a) "
  shows "L - Tmax  exp (-(a * t)) * (L - T)"
    and "L - exp (-(a * t)) * (L - T)  Tmax"
    and "Tmin  L - exp (-(a * t)) * (L - T)"
proof-
  have "0  t  t  - (ln ((L - Tmax) / (L - T)) / a)"
    using thyps by auto
  hence "ln ((L - Tmax) / (L - T))  - a * t  - a * t  0"
    using assms(1) divide_le_cancel by fastforce
  also have "(L - Tmax) / (L - T) > 0"
    using Thyps by auto
  ultimately have "(L - Tmax) / (L - T)  exp (-a * t)  exp (-a * t)  1"
    using exp_ln exp_le_one_iff by (metis exp_less_cancel_iff not_less)
  moreover have "L - T > 0"
    using Thyps by auto
  ultimately have obs: "(L - Tmax)  exp (-a * t) * (L - T)  exp (-a * t) * (L - T)  (L - T)"
    by (simp add: pos_divide_le_eq)
  thus "(L - Tmax)  exp (-(a * t)) * (L - T)"
    by auto
  thus "L - exp (-(a * t)) * (L - T)  Tmax"
    by auto
  show "Tmin  L - exp (-(a * t)) * (L - T)"
    using Thyps and obs by auto
qed

lemmas fbox_temp_dyn = local_flow.fbox_g_ode_subset[OF local_flow_temp]

lemma thermostat:
  assumes "a > 0" and "0 < Tmin" and "Tmax < L"
  shows "(λs. Tmin  s$1  s$1  Tmax  s$4 = 0) 
  |LOOP
    ― ‹control›
    ((2 ::= (λs. 0));(3 ::= (λs. s$1));
    (IF (λs. s$4 = 0  s$3  Tmin + 1) THEN (4 ::= (λs.1)) ELSE
    (IF (λs. s$4 = 1  s$3  Tmax - 1) THEN (4 ::= (λs.0)) ELSE skip));
    ― ‹dynamics›
    (IF (λs. s$4 = 0) THEN (x´= f a 0 & (λs. s$2  - (ln (Tmin/s$3))/a))
    ELSE (x´= f a L & (λs. s$2  - (ln ((L-Tmax)/(L-s$3)))/a))) )
  INV (λs. Tmin s$1  s$1  Tmax  (s$4 = 0  s$4 = 1))]
  (λs. Tmin  s$1  s$1  Tmax)"
  apply(rule fbox_loopI, simp_all add: fbox_temp_dyn[OF assms(1)] le_fun_def)
  using temp_dyn_up_real_arith[OF assms(1) _ _ assms(3), of Tmin]
    and temp_dyn_down_real_arith[OF assms(1,2), of _ Tmax] by auto

no_notation temp_vec_field ("f")
        and temp_flow ("φ")

subsubsection ‹ Tank ›

text ‹ A controller turns a water pump on and off to keep the level of water @{text "h"} in a tank
within an acceptable range @{text "hmin ≤ h ≤ hmax"}. Just like in the previous example, after 
each intervention, the controller registers the current level of water and resets its chronometer,
then it changes the status of the water pump accordingly. The level of water grows linearly 
@{text "h' = k"} at a rate of @{text "k = ci-co"} if the pump is on, and at a rate of 
@{text "k = -co"} if the pump is off. We use @{term "1::4"} to denote the tank's level of water,
@{term "2::4"} is time as measured by the controller's chronometer, @{term "3::4"} is the
level of water measured by the chronometer, and @{term "4::4"} states whether the pump is on
(@{text "s$4 = 1"}) or off (@{text "s$4 = 0"}). We prove that the controller keeps the level of
water between @{text "hmin"} and @{text "hmax"}. ›

abbreviation tank_vec_field :: "real  real^4  real^4" ("f")
  where "f k s  (χ i. if i = 2 then 1 else (if i = 1 then k else 0))"

abbreviation tank_flow :: "real  real  real^4  real^4" ("φ")
  where "φ k τ s  (χ i. if i = 1 then k * τ + s$1 else 
  (if i = 2 then τ + s$2 else s$i))"

abbreviation tank_guard :: "real  real  real^4  bool" ("G")
  where "G Hm k s  s$2  (Hm - s$3)/k"

abbreviation tank_loop_inv :: "real  real  real^4  bool" ("I")
  where "I hmin hmax s  hmin  s$1  s$1  hmax  (s$4 = 0  s$4 = 1)"

abbreviation tank_diff_inv :: "real  real  real  real^4  bool" ("dI")
  where "dI hmin hmax k s  s$1 = k * s$2 + s$3  0  s$2  
    hmin  s$3  s$3  hmax  (s$4 =0  s$4 = 1)"

lemma local_flow_tank: "local_flow (f k) UNIV UNIV (φ k)"
  apply (unfold_locales, unfold local_lipschitz_def lipschitz_on_def, simp_all, clarsimp)
  apply(rule_tac x="1/2" in exI, clarsimp, rule_tac x=1 in exI)
  apply(simp add: dist_norm norm_vec_def L2_set_def, unfold UNIV_4)
  by (auto intro!: poly_derivatives simp: vec_eq_iff)

lemma tank_arith:
  assumes "0  (τ::real)" and "0 < co" and "co < ci"
  shows "τ{0..τ}. τ  - ((hmin - y) / co)   hmin  y - co * τ"
    and "τ{0..τ}. τ  (hmax - y) / (ci - co)   (ci - co) * τ + y  hmax"
    and "hmin  y  hmin  (ci - co) * τ + y"
    and "y  hmax  y - co * τ  hmax"
  apply(simp_all add: field_simps le_divide_eq assms)
  using assms apply (meson add_mono less_eq_real_def mult_left_mono)
  using assms by (meson add_increasing2 less_eq_real_def mult_nonneg_nonneg) 

lemma tank_flow:
  assumes "0 < co" and "co < ci"
  shows "I hmin hmax 
  |LOOP 
    ― ‹control›
    ((2 ::=(λs.0));(3 ::=(λs. s$1));
    (IF (λs. s$4 = 0  s$3  hmin + 1) THEN (4 ::= (λs.1)) ELSE 
    (IF (λs. s$4 = 1  s$3  hmax - 1) THEN (4 ::= (λs.0)) ELSE skip));
    ― ‹dynamics›
    (IF (λs. s$4 = 0) THEN (x´= f (ci-co) & (G hmax (ci-co))) 
     ELSE (x´= f (-co) & (G hmin (-co)))) ) INV I hmin hmax]  
  I hmin hmax"
  apply(rule fbox_loopI, simp_all add: le_fun_def)
  apply(clarsimp simp: le_fun_def local_flow.fbox_g_ode_subset[OF local_flow_tank])
  using assms tank_arith[OF _ assms] by auto

no_notation tank_vec_field ("f")
        and tank_flow ("φ")
        and tank_loop_inv ("I")
        and tank_diff_inv ("dI")
        and tank_guard ("G")

end