Skip to content

Stack optimizer constraints so SLSQP differentiates them once per point #42

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Ickaser merged 3 commits into from
Aug 6, 2026
Merged

Stack optimizer constraints so SLSQP differentiates them once per point #42

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
84630c3
Stack optimizer constraints so SLSQP differentiates them once per point
bernalde Aug 5, 2026
10e5e6e
Rename objective function and its Jacobian in opt_Pch_Tsh
Ickaser Aug 6, 2026
acd560c
Rename objective function and its Jacobian in opt_Tsh
Ickaser Aug 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 18 additions & 10 deletions lyopronto/opt_Pch.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
29 changes: 19 additions & 10 deletions lyopronto/opt_Pch_Tsh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,21 +52,30 @@ 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 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
# 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(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
Expand Down
30 changes: 19 additions & 11 deletions lyopronto/opt_Tsh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,22 +59,30 @@ 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 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
# 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(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
Expand Down
Loading