Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py: 8%
104 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Byrd-Omojokun Trust-Region SQP method."""
3from scipy.sparse import eye as speye
4from .projections import projections
5from .qp_subproblem import modified_dogleg, projected_cg, box_intersections
6import numpy as np
7from numpy.linalg import norm
9__all__ = ['equality_constrained_sqp']
12def default_scaling(x):
13 n, = np.shape(x)
14 return speye(n)
17def equality_constrained_sqp(fun_and_constr, grad_and_jac, lagr_hess,
18 x0, fun0, grad0, constr0,
19 jac0, stop_criteria,
20 state,
21 initial_penalty,
22 initial_trust_radius,
23 factorization_method,
24 trust_lb=None,
25 trust_ub=None,
26 scaling=default_scaling):
27 """Solve nonlinear equality-constrained problem using trust-region SQP.
29 Solve optimization problem:
31 minimize fun(x)
32 subject to: constr(x) = 0
34 using Byrd-Omojokun Trust-Region SQP method described in [1]_. Several
35 implementation details are based on [2]_ and [3]_, p. 549.
37 References
38 ----------
39 .. [1] Lalee, Marucha, Jorge Nocedal, and Todd Plantenga. "On the
40 implementation of an algorithm for large-scale equality
41 constrained optimization." SIAM Journal on
42 Optimization 8.3 (1998): 682-706.
43 .. [2] Byrd, Richard H., Mary E. Hribar, and Jorge Nocedal.
44 "An interior point algorithm for large-scale nonlinear
45 programming." SIAM Journal on Optimization 9.4 (1999): 877-900.
46 .. [3] Nocedal, Jorge, and Stephen J. Wright. "Numerical optimization"
47 Second Edition (2006).
48 """
49 PENALTY_FACTOR = 0.3 # Rho from formula (3.51), reference [2]_, p.891.
50 LARGE_REDUCTION_RATIO = 0.9
51 INTERMEDIARY_REDUCTION_RATIO = 0.3
52 SUFFICIENT_REDUCTION_RATIO = 1e-8 # Eta from reference [2]_, p.892.
53 TRUST_ENLARGEMENT_FACTOR_L = 7.0
54 TRUST_ENLARGEMENT_FACTOR_S = 2.0
55 MAX_TRUST_REDUCTION = 0.5
56 MIN_TRUST_REDUCTION = 0.1
57 SOC_THRESHOLD = 0.1
58 TR_FACTOR = 0.8 # Zeta from formula (3.21), reference [2]_, p.885.
59 BOX_FACTOR = 0.5
61 n, = np.shape(x0) # Number of parameters
63 # Set default lower and upper bounds.
64 if trust_lb is None:
65 trust_lb = np.full(n, -np.inf)
66 if trust_ub is None:
67 trust_ub = np.full(n, np.inf)
69 # Initial values
70 x = np.copy(x0)
71 trust_radius = initial_trust_radius
72 penalty = initial_penalty
73 # Compute Values
74 f = fun0
75 c = grad0
76 b = constr0
77 A = jac0
78 S = scaling(x)
79 # Get projections
80 Z, LS, Y = projections(A, factorization_method)
81 # Compute least-square lagrange multipliers
82 v = -LS.dot(c)
83 # Compute Hessian
84 H = lagr_hess(x, v)
86 # Update state parameters
87 optimality = norm(c + A.T.dot(v), np.inf)
88 constr_violation = norm(b, np.inf) if len(b) > 0 else 0
89 cg_info = {'niter': 0, 'stop_cond': 0,
90 'hits_boundary': False}
92 last_iteration_failed = False
93 while not stop_criteria(state, x, last_iteration_failed,
94 optimality, constr_violation,
95 trust_radius, penalty, cg_info):
96 # Normal Step - `dn`
97 # minimize 1/2*||A dn + b||^2
98 # subject to:
99 # ||dn|| <= TR_FACTOR * trust_radius
100 # BOX_FACTOR * lb <= dn <= BOX_FACTOR * ub.
101 dn = modified_dogleg(A, Y, b,
102 TR_FACTOR*trust_radius,
103 BOX_FACTOR*trust_lb,
104 BOX_FACTOR*trust_ub)
106 # Tangential Step - `dt`
107 # Solve the QP problem:
108 # minimize 1/2 dt.T H dt + dt.T (H dn + c)
109 # subject to:
110 # A dt = 0
111 # ||dt|| <= sqrt(trust_radius**2 - ||dn||**2)
112 # lb - dn <= dt <= ub - dn
113 c_t = H.dot(dn) + c
114 b_t = np.zeros_like(b)
115 trust_radius_t = np.sqrt(trust_radius**2 - np.linalg.norm(dn)**2)
116 lb_t = trust_lb - dn
117 ub_t = trust_ub - dn
118 dt, cg_info = projected_cg(H, c_t, Z, Y, b_t,
119 trust_radius_t,
120 lb_t, ub_t)
122 # Compute update (normal + tangential steps).
123 d = dn + dt
125 # Compute second order model: 1/2 d H d + c.T d + f.
126 quadratic_model = 1/2*(H.dot(d)).dot(d) + c.T.dot(d)
127 # Compute linearized constraint: l = A d + b.
128 linearized_constr = A.dot(d)+b
129 # Compute new penalty parameter according to formula (3.52),
130 # reference [2]_, p.891.
131 vpred = norm(b) - norm(linearized_constr)
132 # Guarantee `vpred` always positive,
133 # regardless of roundoff errors.
134 vpred = max(1e-16, vpred)
135 previous_penalty = penalty
136 if quadratic_model > 0:
137 new_penalty = quadratic_model / ((1-PENALTY_FACTOR)*vpred)
138 penalty = max(penalty, new_penalty)
139 # Compute predicted reduction according to formula (3.52),
140 # reference [2]_, p.891.
141 predicted_reduction = -quadratic_model + penalty*vpred
143 # Compute merit function at current point
144 merit_function = f + penalty*norm(b)
145 # Evaluate function and constraints at trial point
146 x_next = x + S.dot(d)
147 f_next, b_next = fun_and_constr(x_next)
148 # Compute merit function at trial point
149 merit_function_next = f_next + penalty*norm(b_next)
150 # Compute actual reduction according to formula (3.54),
151 # reference [2]_, p.892.
152 actual_reduction = merit_function - merit_function_next
153 # Compute reduction ratio
154 reduction_ratio = actual_reduction / predicted_reduction
156 # Second order correction (SOC), reference [2]_, p.892.
157 if reduction_ratio < SUFFICIENT_REDUCTION_RATIO and \
158 norm(dn) <= SOC_THRESHOLD * norm(dt):
159 # Compute second order correction
160 y = -Y.dot(b_next)
161 # Make sure increment is inside box constraints
162 _, t, intersect = box_intersections(d, y, trust_lb, trust_ub)
163 # Compute tentative point
164 x_soc = x + S.dot(d + t*y)
165 f_soc, b_soc = fun_and_constr(x_soc)
166 # Recompute actual reduction
167 merit_function_soc = f_soc + penalty*norm(b_soc)
168 actual_reduction_soc = merit_function - merit_function_soc
169 # Recompute reduction ratio
170 reduction_ratio_soc = actual_reduction_soc / predicted_reduction
171 if intersect and reduction_ratio_soc >= SUFFICIENT_REDUCTION_RATIO:
172 x_next = x_soc
173 f_next = f_soc
174 b_next = b_soc
175 reduction_ratio = reduction_ratio_soc
177 # Readjust trust region step, formula (3.55), reference [2]_, p.892.
178 if reduction_ratio >= LARGE_REDUCTION_RATIO:
179 trust_radius = max(TRUST_ENLARGEMENT_FACTOR_L * norm(d),
180 trust_radius)
181 elif reduction_ratio >= INTERMEDIARY_REDUCTION_RATIO:
182 trust_radius = max(TRUST_ENLARGEMENT_FACTOR_S * norm(d),
183 trust_radius)
184 # Reduce trust region step, according to reference [3]_, p.696.
185 elif reduction_ratio < SUFFICIENT_REDUCTION_RATIO:
186 trust_reduction = ((1-SUFFICIENT_REDUCTION_RATIO) /
187 (1-reduction_ratio))
188 new_trust_radius = trust_reduction * norm(d)
189 if new_trust_radius >= MAX_TRUST_REDUCTION * trust_radius:
190 trust_radius *= MAX_TRUST_REDUCTION
191 elif new_trust_radius >= MIN_TRUST_REDUCTION * trust_radius:
192 trust_radius = new_trust_radius
193 else:
194 trust_radius *= MIN_TRUST_REDUCTION
196 # Update iteration
197 if reduction_ratio >= SUFFICIENT_REDUCTION_RATIO:
198 x = x_next
199 f, b = f_next, b_next
200 c, A = grad_and_jac(x)
201 S = scaling(x)
202 # Get projections
203 Z, LS, Y = projections(A, factorization_method)
204 # Compute least-square lagrange multipliers
205 v = -LS.dot(c)
206 # Compute Hessian
207 H = lagr_hess(x, v)
208 # Set Flag
209 last_iteration_failed = False
210 # Otimality values
211 optimality = norm(c + A.T.dot(v), np.inf)
212 constr_violation = norm(b, np.inf) if len(b) > 0 else 0
213 else:
214 penalty = previous_penalty
215 last_iteration_failed = True
217 return x, state