feat: truss design case study

Author: ramonfmir

CvxLean/Command/Equivalence.lean

@@ -111,6 +111,16 @@ def evalEquivalenceAux (probIdStx eqvIdStx : TSyntax `ident) (xs : Array (Syntax
       simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.map_domain
       simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.add_constraint
       simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_objFun
+      simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_objFun_1
+      simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_objFun_2
+      simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_objFun_3
+      simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_objFun_4
+      simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_objFun_5
+      simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_objFun_6
+      simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_objFun_7
+      simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_objFun_8
+      simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_objFun_9
+      simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_objFun_10
       simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_constraints
       simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_constraint_1
       simpThms ← simpThms.addDeclToUnfold ``Minimization.Equivalence.rewrite_constraint_2
@@ -144,6 +154,7 @@ def evalEquivalenceAux (probIdStx eqvIdStx : TSyntax `ident) (xs : Array (Syntax
       let eqvNonPropArgs ← eqvArgs.filterM fun arg => do
         return !(← inferType (← inferType arg)).isProp
       let psi ← mkLambdaFVars eqvNonPropArgs res.expr
+      trace[Meta.debug] "psi: {psi}"
 
       try
         let psiF ← realToFloat psi

CvxLean/Examples/FittingSphere.lean

@@ -80,10 +80,10 @@ def fittingSphere :=
     subject to
       h₁ : 0 < r
 
-instance : ChangeOfVariables fun (ct : (Fin n → ℝ) × ℝ) => (ct.1, sqrt (ct.2 + ‖ct.1‖ ^ 2)) :=
+instance : ChangeOfVariables fun ((c, t) : (Fin n → ℝ) × ℝ) => (c, sqrt (t + ‖c‖ ^ 2)) :=
   { inv := fun (c, r) => (c, r ^ 2 - ‖c‖ ^ 2),
-    condition := fun (_, t) => 0 ≤ t,
-    property := fun ⟨c, t⟩ h => by simp [sqrt_sq h] }
+    condition := fun (_, r) => 0 ≤ r,
+    property := fun ⟨c, r⟩ h => by simp [sqrt_sq h] }
 
 equivalence' eqv/fittingSphereT (n m : ℕ) (x : Fin m → Fin n → ℝ) : fittingSphere n m x := by
   -- Change of variables.

CvxLean/Examples/HypersonicShapeDesign.lean

@@ -1,4 +1,5 @@
 import CvxLean
+import CvxLean.Command.Util.TimeCmd
 
 noncomputable section
 
@@ -16,50 +17,109 @@ 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
+      h₁ : 10e-6 ≤ Δx
+      h₂ : Δx ≤ 1
+      h₃ : 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) :
+equivalence' eqv₁/hypersonicShapeDesignConvex (a b : ℝ) (ha : 0 ≤ a) (hb₁ : 0 ≤ b) (hb₂ : b < 1) :
     hypersonicShapeDesign a b := by
   pre_dcp
 
+#print hypersonicShapeDesignConvex
+
 @[optimization_param]
 def aₚ : ℝ := 0.05
 
+@[simp high]
 lemma aₚ_nonneg : 0 ≤ aₚ := by
   unfold aₚ; norm_num
 
 @[optimization_param]
 def bₚ : ℝ := 0.65
 
+lemma bₚ_nonneg : 0 ≤ bₚ := by
+  unfold bₚ; norm_num
+
 lemma bₚ_lt_one : bₚ < 1 := by
   unfold bₚ; norm_num
 
 @[simp high]
-lemma one_sub_bₚ_nonpos : 0 ≤ 1 - bₚ := by
+lemma one_sub_bₚ_nonneg : 0 ≤ 1 - bₚ := by
   unfold bₚ; norm_num
 
-solve hypersonicShapeDesignConvex aₚ bₚ aₚ_nonneg bₚ_lt_one
+time_cmd solve hypersonicShapeDesignConvex aₚ bₚ aₚ_nonneg bₚ_nonneg bₚ_lt_one
+
+#print hypersonicShapeDesignConvex.reduced
 
 -- Final width of wedge.
-def width := hypersonicShapeDesignConvex.solution
+def width := eqv₁.backward_map aₚ.float bₚ.float hypersonicShapeDesignConvex.solution
 
-#eval width
+#eval width -- 0.989524
+
+#eval aₚ.float * (1 / width) - (1 - bₚ.float) * Float.sqrt (1 - width ^ 2) ≤ 0
+#eval aₚ.float * (1 / width) - (1 - bₚ.float) * Float.sqrt (1 - width ^ 2) ≤ 0.000001
 
 -- Final height of wedge.
 def height := Float.sqrt (1 - width ^ 2)
 
-#eval height
+#eval height -- 0.144368
 
 -- Final L/D ratio.
 def ldRatio := 1 / (Float.sqrt ((1 / width ^ 2) - 1))
 
-#eval ldRatio
+#eval ldRatio -- 6.854156
+
+-- While the above is good enough, we simplify the problem further by performing a change of
+-- variables and simplifying appropriately.
+
+equivalence' eqv₂/hypersonicShapeDesignSimpler (a b : ℝ) (ha : 0 ≤ a) (hb₁ : 0 ≤ b)
+    (hb₂ : b < 1) : hypersonicShapeDesignConvex a b ha hb₁ hb₂ := by
+  change_of_variables (z) (Δx ↦ sqrt z)
+  conv_constr h₁ =>
+    rewrite [le_sqrt' (by norm_num)]; norm_num
+  conv_constr h₂ =>
+    rewrite [sqrt_le_iff]; norm_num
+  rw_constr h₃ =>
+    have hz : 0 ≤ z := by arith
+    have h_one_sub_z : 0 ≤ 1 - z := by arith
+    rewrite [rpow_two (sqrt a), sq_sqrt ha, rpow_two (sqrt z), sq_sqrt hz]
+    rewrite [div_le_iff (by arith)]
+    have hlhs : 0 ≤ a / sqrt z := div_nonneg ha (sqrt_nonneg _)
+    have hrhs : 0 ≤ sqrt (1 - z) * (1 - b) := mul_nonneg (sqrt_nonneg _) (by arith)
+    rewrite [← pow_two_le_pow_two hlhs hrhs]
+    rewrite [div_rpow ha (sqrt_nonneg _), rpow_two (sqrt z), sq_sqrt hz]
+    rewrite [mul_rpow (sqrt_nonneg _) (by arith), rpow_two (sqrt (1 - z)), sq_sqrt h_one_sub_z]
+    rewrite [← mul_one_div, ← inv_eq_one_div, mul_comm (1 - z) _]
+    rfl
+  rename_constrs [h₁, h₂, h₃]
+  rw_obj =>
+    rewrite [rpow_neg (sqrt_nonneg _), rpow_two (sqrt z), sq_sqrt (by arith)]
+    rfl
+
+#print hypersonicShapeDesignSimpler
+
+time_cmd solve hypersonicShapeDesignSimpler aₚ bₚ aₚ_nonneg bₚ_nonneg bₚ_lt_one
+
+#print hypersonicShapeDesignSimpler.reduced
+
+-- Final width of wedge.
+def width' :=
+  eqv₁.backward_map aₚ.float bₚ.float <|
+    eqv₂.backward_map aₚ.float bₚ.float hypersonicShapeDesignSimpler.solution
+
+#eval width' -- 0.989524
+
+#eval aₚ.float * (1 / width') - (1 - bₚ.float) * Float.sqrt (1 - width' ^ 2) ≤ 0
+
+-- Final height of wedge.
+def height' := Float.sqrt (1 - width' ^ 2)
+
+#eval height' -- 0.144371
+
+-- Final L/D ratio.
+def ldRatio' := 1 / (Float.sqrt ((1 / width' ^ 2) - 1))
+
+#eval ldRatio' -- 6.854031
 
 end HypersonicShapeDesign