From 84630c33a403a5d5a3642bca64f38f6e79cd73d1 Mon Sep 17 00:00:00 2001 From: "David E. Bernal Neira" Date: Wed, 5 Aug 2026 15:54:08 -0400 Subject: [PATCH 1/3] Stack optimizer constraints so SLSQP differentiates them once per point Each component of the constraint system was its own entry in the cons tuple, so Eq_Constraints computed all four residuals on every call while each lambda discarded three, and SLSQP finite-differenced each entry separately -- a full 8-point sweep per entry on a 7-variable problem instead of one sweep for the system. dry() solves one such problem per timestep, so the overhead landed on every simulated hour. Stack the equality residuals into one vector-valued constraint and the inequality residuals into another, and supply the exact gradient of the linear objective so SLSQP does not finite-difference that either. Measured on the standard fixtures from tests/conftest.py, median of 3: current stacked + objective gradient opt_Tsh 11.94 s 4.61 s 3.48 s (3.4x) opt_Pch 1.26 s 0.48 s 0.37 s (3.4x) opt_Pch_Tsh 13.42 s 5.55 s 4.14 s (3.2x) Output trajectories are bitwise identical to the previous code on all three cases: same shape, np.array_equal true across all seven columns. The objective is linear, so its finite-difference gradient was already exact; supplying it analytically only removes the evaluations. Full test suite: 144 passed, 1 skipped, in 191 s against 608 s on the parent commit. --- lyopronto/opt_Pch.py | 28 ++++++++++++++++++---------- lyopronto/opt_Pch_Tsh.py | 27 ++++++++++++++++++--------- lyopronto/opt_Tsh.py | 28 ++++++++++++++++++---------- 3 files changed, 54 insertions(+), 29 deletions(-) diff --git a/lyopronto/opt_Pch.py b/lyopronto/opt_Pch.py index 6a9be783..c272670c 100644 --- a/lyopronto/opt_Pch.py +++ b/lyopronto/opt_Pch.py @@ -63,6 +63,10 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial): # Objective function to be minimized to maximize sublimation rate def objfun(x): return (x[0]-x[4]) + # Exact gradient of the linear objective, so SLSQP does not + # finite-difference it at every point. + def objfun_jac(x): + return np.array([1.0,0.0,0.0,0.0,-1.0,0.0,0.0]) # Quantities solved for: x = [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] x0 = np.array([P0,0.0,Tb0,Tsh0,P0*1.1,Ts0,3.0e-4]) # Initial values failures = 0 @@ -71,19 +75,23 @@ def objfun(x): Rp = functions.Rp_FUN(Lck,product['R0'],product['A1'],product['A2']) # Product resistance [cm^2-hr-Torr/g] - # Constraints - cons = ({'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[0]}, # sublimation front pressure [Torr] - {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[1]}, # sublimation rate [kg/hr] - {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[2]}, # vial heat transfer balance - {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[3]}, # shelf temperature [degC] - {'type':'eq','fun':lambda x: x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0])}, # vial heat transfer coefficient [cal/s/K/cm^2] - {'type':'eq','fun':lambda x: x[3]-Tsh}, # shelf temperature fixed [degC] - {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[0]}, # equipment capability inequlity - {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[1]}) # maximum product temperature inequality + # Stack the equality constraints into one vector-valued constraint so + # SLSQP evaluates and differentiates the whole system once per point + # rather than once per component: sublimation front pressure [Torr], + # sublimation rate [kg/hr], vial heat transfer balance, shelf + # temperature [degC], vial heat transfer coefficient [cal/s/K/cm^2], and fixed shelf temperature [degC] + def eq_sys(x, Tsh=Tsh): + return np.array(functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp) + + (x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0]), x[3]-Tsh)) + # Inequality constraints: equipment capability and maximum product temperature + def ineq_sys(x): + return np.array(functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)) + cons = ({'type':'eq','fun':eq_sys}, + {'type':'ineq','fun':ineq_sys}) # Bounds for the unknowns bnds = ((Pchamber['min'],Pchamber.get('max', None)),(0,None),(None,None),(None,None),(0,None),(None,None),(0,None)) # Minimize the objective function i.e. maximize the sublimation rate - res = sp.minimize(objfun,x0,bounds = bnds, constraints = cons) + res = sp.minimize(objfun,x0,jac = objfun_jac,bounds = bnds, constraints = cons) [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] = res['x'] # Results [Torr], [kg/hr], [degC], [degC], [Torr], [degC], [cal/s/K/cm^2] # # Use the results as a guess for the next iteration # TODO: decide on appropriate error handling for unsuccessful iterations diff --git a/lyopronto/opt_Pch_Tsh.py b/lyopronto/opt_Pch_Tsh.py index 214edf81..e836a557 100644 --- a/lyopronto/opt_Pch_Tsh.py +++ b/lyopronto/opt_Pch_Tsh.py @@ -54,19 +54,28 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial): # Quantities solved for: x = [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] def fun(x): return x[0]-x[4] # Objective function to be minimized to maximize sublimation rate + # Exact gradient of the linear objective, so SLSQP does not + # finite-difference it at every point. + def fun_jac(x): + return np.array([1.0,0.0,0.0,0.0,-1.0,0.0,0.0]) x0 = [P0,0.0,T0,T0,P0,T0,3.0e-4] # Initial values - # Constraints - cons = ({'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[0]}, # sublimation front pressure [Torr] - {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[1]}, # sublimation rate [kg/hr] - {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[2]}, # vial heat transfer balance - {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[3]}, # shelf temperature [degC] - {'type':'eq','fun':lambda x: x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0])}, # vial heat transfer coefficient [cal/s/K/cm^2] - {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[0]}, # equipment capability inequlity - {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[1]}) # maximum product temperature inequality + # Stack the equality constraints into one vector-valued constraint so + # SLSQP evaluates and differentiates the whole system once per point + # rather than once per component: sublimation front pressure [Torr], + # sublimation rate [kg/hr], vial heat transfer balance, shelf + # temperature [degC], vial heat transfer coefficient [cal/s/K/cm^2] + def eq_sys(x): + return np.array(functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp) + + (x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0]),)) + # Inequality constraints: equipment capability and maximum product temperature + def ineq_sys(x): + return np.array(functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)) + cons = ({'type':'eq','fun':eq_sys}, + {'type':'ineq','fun':ineq_sys}) # Bounds for the unknowns bnds = ((Pchamber['min'],Pchamber.get('max', None)),(None,None),(None,None),(Tshelf['min'],Tshelf['max']),(None,None),(None,None),(None,None)) # Minimize the objective function i.e. maximize the sublimation rate - res = sp.minimize(fun,x0,bounds = bnds, constraints = cons) + res = sp.minimize(fun,x0,jac = fun_jac,bounds = bnds, constraints = cons) [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] = res['x'] # Results [Torr], [kg/hr], [degC], [degC], [Torr], [degC], [cal/s/K/cm^2] # Sublimated ice length diff --git a/lyopronto/opt_Tsh.py b/lyopronto/opt_Tsh.py index 20f6b3e9..9bd959bb 100644 --- a/lyopronto/opt_Tsh.py +++ b/lyopronto/opt_Tsh.py @@ -61,20 +61,28 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial): # Quantities solved for: x = [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] def fun(x): return (x[0]-x[4]) # Objective function to be minimized to maximize sublimation rate + # Exact gradient of the linear objective, so SLSQP does not + # finite-difference it at every point. + def fun_jac(x): + return np.array([1.0,0.0,0.0,0.0,-1.0,0.0,0.0]) x0 = [Pch,0.0,T0,T0,Pch,T0,3.0e-4] # Initial values - # Constraints - cons = ({'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[0]}, # sublimation front pressure [Torr] - {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[1]}, # sublimation rate [kg/hr] - {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[2]}, # vial heat transfer balance - {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[3]}, # shelf temperature [degC] - {'type':'eq','fun':lambda x: x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0])}, # vial heat transfer coefficient [cal/s/K/cm^2] - {'type':'eq','fun':lambda x: x[0]-Pch}, # chamber pressure fixed [Torr] - {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[0]}, # equipment capability inequlity - {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[1]}) # maximum product temperature inequality + # Stack the equality constraints into one vector-valued constraint so + # SLSQP evaluates and differentiates the whole system once per point + # rather than once per component: sublimation front pressure [Torr], + # sublimation rate [kg/hr], vial heat transfer balance, shelf + # temperature [degC], vial heat transfer coefficient [cal/s/K/cm^2], and fixed chamber pressure [Torr] + def eq_sys(x, Pch=Pch): + return np.array(functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp) + + (x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0]), x[0]-Pch)) + # Inequality constraints: equipment capability and maximum product temperature + def ineq_sys(x): + return np.array(functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)) + cons = ({'type':'eq','fun':eq_sys}, + {'type':'ineq','fun':ineq_sys}) # Bounds for the unknowns bnds = ((None,None),(None,None),(None,None),(Tshelf['min'],Tshelf['max']),(None,None),(None,None),(None,None)) # Minimize the objective function i.e. maximize the sublimation rate - res = sp.minimize(fun,x0,bounds = bnds, constraints = cons) + res = sp.minimize(fun,x0,jac = fun_jac,bounds = bnds, constraints = cons) [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] = res['x'] # Results [Torr], [kg/hr], [degC], [degC], [Torr], [degC], [cal/s/K/cm^2] # Sublimated ice length From 10e5e6e7456848ce551895e08e0b1e2c3975df41 Mon Sep 17 00:00:00 2001 From: Isaac Wheeler <47340776+Ickaser@users.noreply.github.com> Date: Wed, 5 Aug 2026 20:56:14 -0400 Subject: [PATCH 2/3] Rename objective function and its Jacobian in opt_Pch_Tsh --- lyopronto/opt_Pch_Tsh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lyopronto/opt_Pch_Tsh.py b/lyopronto/opt_Pch_Tsh.py index e836a557..23b4ff83 100644 --- a/lyopronto/opt_Pch_Tsh.py +++ b/lyopronto/opt_Pch_Tsh.py @@ -52,11 +52,11 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial): Rp = functions.Rp_FUN(Lck,product['R0'],product['A1'],product['A2']) # Product resistance [cm^2-hr-Torr/g] # Quantities solved for: x = [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] - def fun(x): + def objfun(x): return x[0]-x[4] # Objective function to be minimized to maximize sublimation rate # Exact gradient of the linear objective, so SLSQP does not # finite-difference it at every point. - def fun_jac(x): + def objfun_jac(x): return np.array([1.0,0.0,0.0,0.0,-1.0,0.0,0.0]) x0 = [P0,0.0,T0,T0,P0,T0,3.0e-4] # Initial values # Stack the equality constraints into one vector-valued constraint so @@ -75,7 +75,7 @@ def ineq_sys(x): # Bounds for the unknowns bnds = ((Pchamber['min'],Pchamber.get('max', None)),(None,None),(None,None),(Tshelf['min'],Tshelf['max']),(None,None),(None,None),(None,None)) # Minimize the objective function i.e. maximize the sublimation rate - res = sp.minimize(fun,x0,jac = fun_jac,bounds = bnds, constraints = cons) + res = sp.minimize(objfun,x0,jac = objfun_jac,bounds = bnds, constraints = cons) [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] = res['x'] # Results [Torr], [kg/hr], [degC], [degC], [Torr], [degC], [cal/s/K/cm^2] # Sublimated ice length From acd560c0b5c6eee888cdae9f41dddfd31666336c Mon Sep 17 00:00:00 2001 From: Isaac Wheeler <47340776+Ickaser@users.noreply.github.com> Date: Wed, 5 Aug 2026 20:57:27 -0400 Subject: [PATCH 3/3] Rename objective function and its Jacobian in opt_Tsh --- lyopronto/opt_Tsh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lyopronto/opt_Tsh.py b/lyopronto/opt_Tsh.py index 9bd959bb..9a88f71c 100644 --- a/lyopronto/opt_Tsh.py +++ b/lyopronto/opt_Tsh.py @@ -59,11 +59,11 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial): Rp = functions.Rp_FUN(Lck,product['R0'],product['A1'],product['A2']) # Product resistance [cm^2-hr-Torr/g] # Quantities solved for: x = [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] - def fun(x): + def objfun(x): return (x[0]-x[4]) # Objective function to be minimized to maximize sublimation rate # Exact gradient of the linear objective, so SLSQP does not # finite-difference it at every point. - def fun_jac(x): + def objfun_jac(x): return np.array([1.0,0.0,0.0,0.0,-1.0,0.0,0.0]) x0 = [Pch,0.0,T0,T0,Pch,T0,3.0e-4] # Initial values # Stack the equality constraints into one vector-valued constraint so @@ -82,7 +82,7 @@ def ineq_sys(x): # Bounds for the unknowns bnds = ((None,None),(None,None),(None,None),(Tshelf['min'],Tshelf['max']),(None,None),(None,None),(None,None)) # Minimize the objective function i.e. maximize the sublimation rate - res = sp.minimize(fun,x0,jac = fun_jac,bounds = bnds, constraints = cons) + res = sp.minimize(objfun,x0,jac = objfun_jac,bounds = bnds, constraints = cons) [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] = res['x'] # Results [Torr], [kg/hr], [degC], [degC], [Torr], [degC], [cal/s/K/cm^2] # Sublimated ice length