feat: hypersonic shape design case study

CvxLean/Examples/CovarianceEstimation.lean

@@ -1,15 +1,8 @@
-import CvxLean.Lib.Math.Data.Vec
-import CvxLean.Lib.Math.CovarianceEstimation
-import CvxLean.Lib.Math.LinearAlgebra.Matrix.PosDef
-import CvxLean.Syntax.Minimization
-import CvxLean.Tactic.DCP.AtomLibrary.All
-import CvxLean.Command.Reduction
-import CvxLean.Command.Solve
+import CvxLean
 
 namespace CovarianceEstimation
 
-open CvxLean Minimization
-open Real BigOperators Matrix
+open CvxLean Minimization Real BigOperators Matrix
 
 noncomputable def problem (n : ℕ) (N : ℕ) (α : ℝ) (y : Fin N → Fin n →  ℝ) :=
   optimization (R : Matrix (Fin n) (Fin n) ℝ)

CvxLean/Examples/FittingSphere.lean

@@ -45,18 +45,18 @@ lemma leastSquares_optimal_eq_mean {n : ℕ} (hn : 0 < n) (a : Fin n → ℝ) (x
     rwa [h_rw_x, h_rw_y, mul_le_mul_left (by positivity), add_le_add_iff_right] at hy
   have hmean := h (mean a)
   simp at hmean
-  have hz := le_antisymm hmean (sq_nonneg _)
-  rwa [sq_eq_zero_iff, sub_eq_zero] at hz
+  have h_sq_eq_zero := le_antisymm hmean (sq_nonneg _)
+  rwa [sq_eq_zero_iff, sub_eq_zero] at h_sq_eq_zero
 
-def Vec.leastSquares {n : ℕ} (a : Fin n → ℝ) :=
+def leastSquaresVec {n : ℕ} (a : Fin n → ℝ) :=
   optimization (x : ℝ)
     minimize (Vec.sum ((a - Vec.const n x) ^ 2) : ℝ)
 
 /-- Same as `leastSquares_optimal_eq_mean` in vector notation. -/
-lemma vec_leastSquares_optimal_eq_mean {n : ℕ} (hn : 0 < n) (a : Fin n → ℝ) (x : ℝ)
-  (h : (Vec.leastSquares a).optimal x) : x = mean a := by
+lemma leastSquaresVec_optimal_eq_mean {n : ℕ} (hn : 0 < n) (a : Fin n → ℝ) (x : ℝ)
+  (h : (leastSquaresVec a).optimal x) : x = mean a := by
   apply leastSquares_optimal_eq_mean hn a
-  simp [Vec.leastSquares, leastSquares, optimal, feasible] at h ⊢
+  simp [leastSquaresVec, leastSquares, optimal, feasible] at h ⊢
   intros y
   simp only [Vec.sum, Pi.pow_apply, Pi.sub_apply, Vec.const, rpow_two] at h
   exact h y
@@ -85,8 +85,6 @@ instance : ChangeOfVariables fun (ct : (Fin n → ℝ) × ℝ) => (ct.1, sqrt (c
     condition := fun (_, t) => 0 ≤ t,
     property := fun ⟨c, t⟩ h => by simp [sqrt_sq h] }
 
-set_option trace.Meta.debug true
-
 equivalence' eqv/fittingSphereT (n m : ℕ) (x : Fin m → Fin n → ℝ) : fittingSphere n m x := by
   -- Change of variables.
   equivalence_step =>
@@ -103,19 +101,21 @@ equivalence' eqv/fittingSphereT (n m : ℕ) (x : Fin m → Fin n → ℝ) : fitt
       (g := fun (ct : (Fin n → ℝ) × ℝ) =>
         Vec.sum (((Vec.norm x) ^ 2 - 2 * (Matrix.mulVec x ct.1) - Vec.const m ct.2) ^ 2))
     . rintro ⟨c, t⟩ h
-      dsimp at h ⊢; simp [Vec.sum, Vec.norm, Vec.const]
-      congr; funext i; congr 1;
-      rw [@norm_sub_sq ℝ (Fin n → ℝ) _ (PiLp.normedAddCommGroup _ _) (PiLp.innerProductSpace _)]
-      rw [sq_sqrt (rpow_two _ ▸ le_of_lt (sqrt_pos.mp <| h))]
+      dsimp at h ⊢; simp [Vec.sum, Vec.norm, Vec.const]; congr; funext i; congr 1;
+      rw [norm_sub_sq (𝕜 := ℝ) (E := Fin n → ℝ), sq_sqrt (rpow_two _ ▸ le_of_lt (sqrt_pos.mp h))]
       simp [mulVec, inner, dotProduct]
   rename_vars [c, t]
 
 #print fittingSphereT
 
-relaxation rel/fittingSphereConvex (n m : ℕ) (x : Fin m → Fin n → ℝ) : fittingSphereT n m x := by
-  relaxation_step =>
-    apply Relaxation.weaken_constraint (cs' := fun _ => True)
-    . rintro ⟨c, t⟩ _; trivial
+-- Next, we proceed to remove the non-convex constraint by arguing that any (non-trivial) point that
+-- minimizes the objective function wihtout the constraint, also satisfies the constraint. We define
+-- the problem directly, bot note that we could also remove the constraint using the `relaxation`
+-- command.
+
+def fittingSphereConvex (n m : ℕ) (x : Fin m → Fin n → ℝ) :=
+  optimization (c : Fin n → ℝ) (t : ℝ)
+    minimize (Vec.sum ((Vec.norm x ^ 2 - 2 * mulVec x c - Vec.const m t) ^ 2) : ℝ)
 
 /-- If the squared error is zero, then `aᵢ = x`. -/
 lemma vec_squared_norm_error_eq_zero_iff {n m : ℕ} (a : Fin m → Fin n → ℝ) (x : Fin n → ℝ) :
@@ -125,54 +125,53 @@ lemma vec_squared_norm_error_eq_zero_iff {n m : ℕ} (a : Fin m → Fin n → 
   constructor
   . intros h i
     have hi := h i (by simp)
-    rw [sq_eq_zero_iff, @norm_eq_zero _ (PiLp.normedAddCommGroup _ _).toNormedAddGroup] at hi
-    rwa [sub_eq_zero] at hi
+    rwa [sq_eq_zero_iff, norm_eq_zero, sub_eq_zero] at hi
   . intros h i _
-    rw [sq_eq_zero_iff, @norm_eq_zero _ (PiLp.normedAddCommGroup _ _).toNormedAddGroup, sub_eq_zero]
+    rw [sq_eq_zero_iff, norm_eq_zero, sub_eq_zero]
     exact h i
 
 /-- This tells us that solving the relaxed problem is sufficient for optimal points if the solution
 is non-trivial. -/
-lemma optimal_relaxed_implies_optimal (hm : 0 < m) (c : Fin n → ℝ) (t : ℝ)
+lemma optimal_convex_implies_optimal_t (hm : 0 < m) (c : Fin n → ℝ) (t : ℝ)
   (h_nontrivial : x ≠ Vec.const m c)
   (h_opt : (fittingSphereConvex n m x).optimal (c, t)) : (fittingSphereT n m x).optimal (c, t) := by
   simp [fittingSphereT, fittingSphereConvex, optimal, feasible] at h_opt ⊢
   constructor
   . let a := Vec.norm x ^ 2 - 2 * mulVec x c
-    have h_ls : optimal (Vec.leastSquares a) t := by
+    have h_ls : optimal (leastSquaresVec a) t := by
       refine ⟨trivial, ?_⟩
       intros y _
-      simp [objFun, Vec.leastSquares]
+      simp [objFun, leastSquaresVec]
       exact h_opt c y
     -- Apply key result about least squares to `a` and `t`.
-    have ht_eq := vec_leastSquares_optimal_eq_mean hm a t h_ls
-    have hc2_eq : ‖c‖ ^ 2 = (1 / m) * ∑ i : Fin m, ‖c‖ ^ 2 := by
+    have h_t_eq := leastSquaresVec_optimal_eq_mean hm a t h_ls
+    have h_c2_eq : ‖c‖ ^ 2 = (1 / m) * ∑ i : Fin m, ‖c‖ ^ 2 := by
       simp [sum_const]
       field_simp; ring
-    have ht : t + ‖c‖ ^ 2 = (1 / m) * ∑ i, ‖(x i) - c‖ ^ 2 := by
-      rw [ht_eq]; dsimp [mean]
-      rw [hc2_eq, mul_sum, mul_sum, mul_sum, ← sum_add_distrib]
+    have h_t_add_c2_eq : t + ‖c‖ ^ 2 = (1 / m) * ∑ i, ‖(x i) - c‖ ^ 2 := by
+      rw [h_t_eq]; dsimp [mean]
+      rw [h_c2_eq, mul_sum, mul_sum, mul_sum, ← sum_add_distrib]
       congr; funext i; rw [← mul_add]
       congr; simp [Vec.norm]
-      rw [@norm_sub_sq ℝ (Fin n → ℝ) _ (PiLp.normedAddCommGroup _ _) (PiLp.innerProductSpace _)]
+      rw [norm_sub_sq (𝕜 := ℝ) (E := Fin n → ℝ)]
       congr
     -- We use the result to establish that `t + ‖c‖ ^ 2` is non-negative.
-    have h_tc2_nonneg : 0 ≤ t + ‖c‖ ^ 2 := by
-      rw [ht]
+    have h_t_add_c2_nonneg : 0 ≤ t + ‖c‖ ^ 2 := by
+      rw [h_t_add_c2_eq]
       apply mul_nonneg (by norm_num)
       apply sum_nonneg
       intros i _
       rw [rpow_two]
       exact sq_nonneg _
-    cases (lt_or_eq_of_le h_tc2_nonneg) with
-    | inl h_tc2_lt_zero =>
+    cases (lt_or_eq_of_le h_t_add_c2_nonneg) with
+    | inl h_t_add_c2_lt_zero =>
         -- If it is positive, we are done.
-        convert h_tc2_lt_zero; simp
-    | inr h_tc2_eq_zero =>
+        convert h_t_add_c2_lt_zero; simp
+    | inr h_t_add_c2_eq_zero =>
         -- Otherwise, it contradicts the non-triviality assumption.
         exfalso
-        rw [ht, zero_eq_mul] at h_tc2_eq_zero
-        rcases h_tc2_eq_zero with (hc | h_sum_eq_zero)
+        rw [h_t_add_c2_eq, zero_eq_mul] at h_t_add_c2_eq_zero
+        rcases h_t_add_c2_eq_zero with (hc | h_sum_eq_zero)
         . simp at hc; linarith
         rw [vec_squared_norm_error_eq_zero_iff] at h_sum_eq_zero
         apply h_nontrivial
@@ -181,6 +180,24 @@ lemma optimal_relaxed_implies_optimal (hm : 0 < m) (c : Fin n → ℝ) (t : ℝ)
   . intros c' x' _
     exact h_opt c' x'
 
+/-- We express the nontriviality condition only in terms of `x` so that it can be checked. -/
+lemma non_triviality_condition (c : Fin n → ℝ) (hx : ∃ i j, x i ≠ x j) : x ≠ Vec.const m c := by
+  intros h
+  conv at hx => congr; ext i; rw [← not_forall]
+  rw [← not_forall] at hx
+  apply hx
+  intros i j
+  rw [congr_fun h i, congr_fun h j]
+  simp [Vec.const]
+
+/-- We show that we have a reduction via the identity map. -/
+def red (hm : 0 < m) (hx : ∃ i j, x i ≠ x j) :
+    (fittingSphereT n m x) ≼ (fittingSphereConvex n m x) :=
+  { psi := id,
+    psi_optimality := fun (c, t) h_opt => by
+      have h_nontrivial := non_triviality_condition n m x c hx
+      exact optimal_convex_implies_optimal_t n m x hm c t h_nontrivial h_opt }
+
 #print fittingSphereConvex
 
 -- We proceed to solve the problem on a concrete example.
@@ -205,16 +222,12 @@ def xₚ : Fin mₚ → Fin nₚ → ℝ := Matrix.transpose <| ![
 
 -- We use the `solve` command on the data above.
 
-set_option maxHeartbeats 1000000
-
 solve fittingSphereConvex nₚ mₚ xₚ
 
 -- Finally, we recover the solution to the original problem.
 
 def sol := eqv.backward_map nₚ mₚ xₚ.float fittingSphereConvex.solution
 
-#print eqv.backward_map
-
 #eval sol -- (![1.664863, 0.031932], 1.159033)
 
 end FittingSphere

CvxLean/Examples/HypersonicShapeDesign.lean

@@ -0,0 +1,66 @@
+import CvxLean
+
+noncomputable section
+
+namespace HypersonicShapeDesign
+
+open CvxLean Minimization Real
+
+-- Height of rectangle.
+variable (a : ℝ)
+
+-- Width of rectangle.
+variable (b : ℝ)
+
+def hypersonicShapeDesign :=
+  optimization (Δx : ℝ)
+    minimize sqrt ((1 / Δx ^ 2) - 1)
+    subject to
+      h1 : 10e-6 ≤ Δx
+      h2 : Δx ≤ 1
+      h3 : a * (1 / Δx) - (1 - b) * sqrt (1 - Δx ^ 2) ≤ 0
+
+set_option trace.Meta.debug true
+
+#check Lean.Expr
+
+equivalence eqv/hypersonicShapeDesignConvex (a b : ℝ) (ha : 0 ≤ a) (hb : b < 1) :
+    hypersonicShapeDesign a b := by
+  pre_dcp
+
+@[optimization_param]
+def aₚ : ℝ := 0.05
+
+lemma aₚ_nonneg : 0 ≤ aₚ := by
+  unfold aₚ; norm_num
+
+@[optimization_param]
+def bₚ : ℝ := 0.65
+
+lemma bₚ_lt_one : bₚ < 1 := by
+  unfold bₚ; norm_num
+
+@[simp high]
+lemma one_sub_bₚ_nonpos : 0 ≤ 1 - bₚ := by
+  unfold bₚ; norm_num
+
+solve hypersonicShapeDesignConvex aₚ bₚ aₚ_nonneg bₚ_lt_one
+
+-- Final width of wedge.
+def width := hypersonicShapeDesignConvex.solution
+
+#eval width
+
+-- Final height of wedge.
+def height := Float.sqrt (1 - width ^ 2)
+
+#eval height
+
+-- Final L/D ratio.
+def ldRatio := 1 / (Float.sqrt ((1 / width ^ 2) - 1))
+
+#eval ldRatio
+
+end HypersonicShapeDesign
+
+end

CvxLean/Examples/VehicleSpeedScheduling.lean

@@ -85,12 +85,12 @@ equivalence' eqv₁/vehSpeedSchedConvex (n : ℕ) (d : Fin n → ℝ)
 
 #check eqv₁.backward_map
 
--- The problem is technically in DCP form. The only issue is that we do not have an atom for the
--- perspective function, so the objective function `Vec.sum (t * Vec.map F (d / t))` cannot be
--- reduced directly.
+-- The problem is technically in DCP form if `F` is DCP convex. The only issue is that we do not
+-- have an atom for the perspective function, so the objective function
+-- `Vec.sum (t * Vec.map F (d / t))` cannot be reduced directly.
 
--- We fix `F` and declare an atom for this particular application of the perspective function.
--- Let `F(s) = a * s^2 + b * s + c` with `0 ≤ a`.
+-- However, if we fix `F`, we can use other atoms. For example, if `F` is quadratic, the problem can
+-- be reduced. Let `F(s) = a * s^2 + b * s + c` with `0 ≤ a`.
 
 equivalence' eqv₂/vehSpeedSchedQuadratic (n : ℕ) (d : Fin n → ℝ)
     (τmin τmax smin smax : Fin n → ℝ) (a b c : ℝ) (h_n_pos : 0 < n) (h_d_pos : StrongLT 0 d)

CvxLean/Examples/fitting_sphere.py

@@ -1,5 +1,6 @@
 import numpy as np
 import cvxpy as cp
+import matplotlib.pyplot as plt
 
 n = 2
 
@@ -30,3 +31,19 @@
 print("t* = ", t.value)
 print("c* = ", c.value)
 print("r* = ", r)
+
+def plot_circle_and_points(center, radius, points):
+    theta = np.linspace(0, 2 * np.pi, 100)
+    x_circle = center[0] + radius * np.cos(theta)
+    y_circle = center[1] + radius * np.sin(theta)
+
+    plt.plot(x_circle, y_circle, label='Circle')
+
+    plt.scatter(points[:, 0], points[:, 1], color='red')
+
+    plt.axis('equal')
+
+    plt.show()
+    plt.savefig('plots/fitting_sphere.png')
+
+plot_circle_and_points(c.value, r, x)

CvxLean/Examples/hypersonic_shape_design.py

@@ -0,0 +1,38 @@
+import cvxpy as cp 
+import numpy as np
+
+a = .05 # height of rectangle
+
+b = .65 # width of rectangle
+
+x = cp.Variable(pos=True)
+
+obj = cp.sqrt(cp.inv_pos(cp.square(x)) - 1)
+
+p = cp.Problem(
+    cp.Minimize(obj), [
+        a * cp.inv_pos(x) - (1 - b) * cp.sqrt(1 - cp.square(x)) <= 0
+    ])
+
+p.solve(qcp=True, verbose=True)
+
+print('QCP Final L/D Ratio = ', 1 / obj.value)
+print('QCP Final width of wedge = ', x.value)
+print('QCP Final height of wedge = ', np.sqrt(1 - x.value ** 2))
+
+x = cp.Variable(pos=True)
+
+obj = cp.power(x, -2) - 1
+
+p = cp.Problem(
+    cp.Minimize(obj), [
+        a * cp.inv_pos(x) - (1 - b) * cp.sqrt(1 - cp.square(x)) <= 0
+    ])
+
+p.solve(verbose=True)
+
+print('DCP Final L/D Ratio = ', 1 / np.sqrt(obj.value))
+print('DCP Final width of wedge = ', x.value)
+print('DCP Final height of wedge = ', np.sqrt(1 - x.value ** 2))
+
+

CvxLean/Examples/vehicle_speed_scheduling.py

@@ -1,5 +1,6 @@
 import numpy as np
 import cvxpy as cp
+import matplotlib.pyplot as plt
 
 n = 10
 
@@ -35,5 +36,14 @@
 
 s = d / t.value
 
+def plot(s):
+    plt.step(np.arange(n), s)
+    plt.xlabel('$i$')
+    plt.ylabel('$s_i$')
+    plt.savefig("plots/vehicle_speed_scheduling.png")
+    plt.show()
+
 print("t* = ", t.value)
 print("s* = ", s)
+
+plot(s)