parameters: (fun_and_constr, grad_and_jac, lagr_hess,
                               x0, fun0, grad0, constr0,
                               jac0, stop_criteria,
                               state,
                               initial_penalty,
                               initial_trust_radius,
                               factorization_method,
                               trust_lb=None,
                               trust_ub=None,
                               scaling=default_scaling)
return_statement: return x, state
docstring: Solve nonlinear equality-constrained problem using trust-region SQP.
  
      Solve optimization problem:
  
          minimize fun(x)
          subject to: constr(x) = 0
  
      using Byrd-Omojokun Trust-Region SQP method described in [1]_. Several
      implementation details are based on [2]_ and [3]_, p. 549.
  
      References
      ----------
      .. [1] Lalee, Marucha, Jorge Nocedal, and Todd Plantenga. "On the
             implementation of an algorithm for large-scale equality
             constrained optimization." SIAM Journal on
             Optimization 8.3 (1998): 682-706.
      .. [2] Byrd, Richard H., Mary E. Hribar, and Jorge Nocedal.
             "An interior point algorithm for large-scale nonlinear
             programming." SIAM Journal on Optimization 9.4 (1999): 877-900.
      .. [3] Nocedal, Jorge, and Stephen J. Wright. "Numerical optimization"
             Second Edition (2006).
docstring_summary: Solve nonlinear equality-constrained problem using trust-region SQP.
docstring_tokens: Solve nonlinear equality - constrained problem using trust - region SQP .
function: def equality_constrained_sqp(fun_and_constr, grad_and_jac, lagr_hess,
                               x0, fun0, grad0, constr0,
                               jac0, stop_criteria,
                               state,
                               initial_penalty,
                               initial_trust_radius,
                               factorization_method,
                               trust_lb=None,
                               trust_ub=None,
                               scaling=default_scaling):
      """Solve nonlinear equality-constrained problem using trust-region SQP.
  
      Solve optimization problem:
  
          minimize fun(x)
          subject to: constr(x) = 0
  
      using Byrd-Omojokun Trust-Region SQP method described in [1]_. Several
      implementation details are based on [2]_ and [3]_, p. 549.
  
      References
      ----------
      .. [1] Lalee, Marucha, Jorge Nocedal, and Todd Plantenga. "On the
             implementation of an algorithm for large-scale equality
             constrained optimization." SIAM Journal on
             Optimization 8.3 (1998): 682-706.
      .. [2] Byrd, Richard H., Mary E. Hribar, and Jorge Nocedal.
             "An interior point algorithm for large-scale nonlinear
             programming." SIAM Journal on Optimization 9.4 (1999): 877-900.
      .. [3] Nocedal, Jorge, and Stephen J. Wright. "Numerical optimization"
             Second Edition (2006).
      """
      PENALTY_FACTOR = 0.3  # Rho from formula (3.51), reference [2]_, p.891.
      LARGE_REDUCTION_RATIO = 0.9
      INTERMEDIARY_REDUCTION_RATIO = 0.3
      SUFFICIENT_REDUCTION_RATIO = 1e-8  # Eta from reference [2]_, p.892.
      TRUST_ENLARGEMENT_FACTOR_L = 7.0
      TRUST_ENLARGEMENT_FACTOR_S = 2.0
      MAX_TRUST_REDUCTION = 0.5
      MIN_TRUST_REDUCTION = 0.1
      SOC_THRESHOLD = 0.1
      TR_FACTOR = 0.8  # Zeta from formula (3.21), reference [2]_, p.885.
      BOX_FACTOR = 0.5
  
      n, = np.shape(x0)  # Number of parameters
  
      # Set default lower and upper bounds.
      if trust_lb is None:
          trust_lb = np.full(n, -np.inf)
      if trust_ub is None:
          trust_ub = np.full(n, np.inf)
  
      # Initial values
      x = np.copy(x0)
      trust_radius = initial_trust_radius
      penalty = initial_penalty
      # Compute Values
      f = fun0
      c = grad0
      b = constr0
      A = jac0
      S = scaling(x)
      # Get projections
      Z, LS, Y = projections(A, factorization_method)
      # Compute least-square lagrange multipliers
      v = -LS.dot(c)
      # Compute Hessian
      H = lagr_hess(x, v)
  
      # Update state parameters
      optimality = norm(c + A.T.dot(v), np.inf)
      constr_violation = norm(b, np.inf) if len(b) > 0 else 0
      cg_info = {'niter': 0, 'stop_cond': 0,
                 'hits_boundary': False}
  
      last_iteration_failed = False
      while not stop_criteria(state, x, last_iteration_failed,
                              optimality, constr_violation,
                              trust_radius, penalty, cg_info):
          # Normal Step - `dn`
          # minimize 1/2*||A dn + b||^2
          # subject to:
          # ||dn|| <= TR_FACTOR * trust_radius
          # BOX_FACTOR * lb <= dn <= BOX_FACTOR * ub.
          dn = modified_dogleg(A, Y, b,
                               TR_FACTOR*trust_radius,
                               BOX_FACTOR*trust_lb,
                               BOX_FACTOR*trust_ub)
  
          # Tangential Step - `dt`
          # Solve the QP problem:
          # minimize 1/2 dt.T H dt + dt.T (H dn + c)
          # subject to:
          # A dt = 0
          # ||dt|| <= sqrt(trust_radius**2 - ||dn||**2)
          # lb - dn <= dt <= ub - dn
          c_t = H.dot(dn) + c
          b_t = np.zeros_like(b)
          trust_radius_t = np.sqrt(trust_radius**2 - np.linalg.norm(dn)**2)
          lb_t = trust_lb - dn
          ub_t = trust_ub - dn
          dt, cg_info = projected_cg(H, c_t, Z, Y, b_t,
                                     trust_radius_t,
                                     lb_t, ub_t)
  
          # Compute update (normal + tangential steps).
          d = dn + dt
  
          # Compute second order model: 1/2 d H d + c.T d + f.
          quadratic_model = 1/2*(H.dot(d)).dot(d) + c.T.dot(d)
          # Compute linearized constraint: l = A d + b.
          linearized_constr = A.dot(d)+b
          # Compute new penalty parameter according to formula (3.52),
          # reference [2]_, p.891.
          vpred = norm(b) - norm(linearized_constr)
          # Guarantee `vpred` always positive,
          # regardless of roundoff errors.
          vpred = max(1e-16, vpred)
          previous_penalty = penalty
          if quadratic_model > 0:
              new_penalty = quadratic_model / ((1-PENALTY_FACTOR)*vpred)
              penalty = max(penalty, new_penalty)
          # Compute predicted reduction according to formula (3.52),
          # reference [2]_, p.891.
          predicted_reduction = -quadratic_model + penalty*vpred
  
          # Compute merit function at current point
          merit_function = f + penalty*norm(b)
          # Evaluate function and constraints at trial point
          x_next = x + S.dot(d)
          f_next, b_next = fun_and_constr(x_next)
          # Compute merit function at trial point
          merit_function_next = f_next + penalty*norm(b_next)
          # Compute actual reduction according to formula (3.54),
          # reference [2]_, p.892.
          actual_reduction = merit_function - merit_function_next
          # Compute reduction ratio
          reduction_ratio = actual_reduction / predicted_reduction
  
          # Second order correction (SOC), reference [2]_, p.892.
          if reduction_ratio < SUFFICIENT_REDUCTION_RATIO and \
             norm(dn) <= SOC_THRESHOLD * norm(dt):
              # Compute second order correction
              y = -Y.dot(b_next)
              # Make sure increment is inside box constraints
              _, t, intersect = box_intersections(d, y, trust_lb, trust_ub)
              # Compute tentative point
              x_soc = x + S.dot(d + t*y)
              f_soc, b_soc = fun_and_constr(x_soc)
              # Recompute actual reduction
              merit_function_soc = f_soc + penalty*norm(b_soc)
              actual_reduction_soc = merit_function - merit_function_soc
              # Recompute reduction ratio
              reduction_ratio_soc = actual_reduction_soc / predicted_reduction
              if intersect and reduction_ratio_soc >= SUFFICIENT_REDUCTION_RATIO:
                  x_next = x_soc
                  f_next = f_soc
                  b_next = b_soc
                  reduction_ratio = reduction_ratio_soc
  
          # Readjust trust region step, formula (3.55), reference [2]_, p.892.
          if reduction_ratio >= LARGE_REDUCTION_RATIO:
              trust_radius = max(TRUST_ENLARGEMENT_FACTOR_L * norm(d),
                                 trust_radius)
          elif reduction_ratio >= INTERMEDIARY_REDUCTION_RATIO:
              trust_radius = max(TRUST_ENLARGEMENT_FACTOR_S * norm(d),
                                 trust_radius)
          # Reduce trust region step, according to reference [3]_, p.696.
          elif reduction_ratio < SUFFICIENT_REDUCTION_RATIO:
                  trust_reduction \
                      = (1-SUFFICIENT_REDUCTION_RATIO)/(1-reduction_ratio)
                  new_trust_radius = trust_reduction * norm(d)
                  if new_trust_radius >= MAX_TRUST_REDUCTION * trust_radius:
                      trust_radius *= MAX_TRUST_REDUCTION
                  elif new_trust_radius >= MIN_TRUST_REDUCTION * trust_radius:
                      trust_radius = new_trust_radius
                  else:
                      trust_radius *= MIN_TRUST_REDUCTION
  
          # Update iteration
          if reduction_ratio >= SUFFICIENT_REDUCTION_RATIO:
              x = x_next
              f, b = f_next, b_next
              c, A = grad_and_jac(x)
              S = scaling(x)
              # Get projections
              Z, LS, Y = projections(A, factorization_method)
              # Compute least-square lagrange multipliers
              v = -LS.dot(c)
              # Compute Hessian
              H = lagr_hess(x, v)
              # Set Flag
              last_iteration_failed = False
              # Otimality values
              optimality = norm(c + A.T.dot(v), np.inf)
              constr_violation = norm(b, np.inf) if len(b) > 0 else 0
          else:
              penalty = previous_penalty
              last_iteration_failed = True
  
      return x, state
function_tokens: def equality_constrained_sqp ( fun_and_constr , grad_and_jac , lagr_hess , x0 , fun0 , grad0 , constr0 , jac0 , stop_criteria , state , initial_penalty , initial_trust_radius , factorization_method , trust_lb = None , trust_ub = None , scaling = default_scaling ) : PENALTY_FACTOR = 0.3 # Rho from formula (3.51), reference [2]_, p.891. LARGE_REDUCTION_RATIO = 0.9 INTERMEDIARY_REDUCTION_RATIO = 0.3 SUFFICIENT_REDUCTION_RATIO = 1e-8 # Eta from reference [2]_, p.892. TRUST_ENLARGEMENT_FACTOR_L = 7.0 TRUST_ENLARGEMENT_FACTOR_S = 2.0 MAX_TRUST_REDUCTION = 0.5 MIN_TRUST_REDUCTION = 0.1 SOC_THRESHOLD = 0.1 TR_FACTOR = 0.8 # Zeta from formula (3.21), reference [2]_, p.885. BOX_FACTOR = 0.5 n , = np . shape ( x0 ) # Number of parameters # Set default lower and upper bounds. if trust_lb is None : trust_lb = np . full ( n , - np . inf ) if trust_ub is None : trust_ub = np . full ( n , np . inf ) # Initial values x = np . copy ( x0 ) trust_radius = initial_trust_radius penalty = initial_penalty # Compute Values f = fun0 c = grad0 b = constr0 A = jac0 S = scaling ( x ) # Get projections Z , LS , Y = projections ( A , factorization_method ) # Compute least-square lagrange multipliers v = - LS . dot ( c ) # Compute Hessian H = lagr_hess ( x , v ) # Update state parameters optimality = norm ( c + A . T . dot ( v ) , np . inf ) constr_violation = norm ( b , np . inf ) if len ( b ) > 0 else 0 cg_info = { 'niter' : 0 , 'stop_cond' : 0 , 'hits_boundary' : False } last_iteration_failed = False while not stop_criteria ( state , x , last_iteration_failed , optimality , constr_violation , trust_radius , penalty , cg_info ) : # Normal Step - `dn` # minimize 1/2*||A dn + b||^2 # subject to: # ||dn|| <= TR_FACTOR * trust_radius # BOX_FACTOR * lb <= dn <= BOX_FACTOR * ub. dn = modified_dogleg ( A , Y , b , TR_FACTOR * trust_radius , BOX_FACTOR * trust_lb , BOX_FACTOR * trust_ub ) # Tangential Step - `dt` # Solve the QP problem: # minimize 1/2 dt.T H dt + dt.T (H dn + c) # subject to: # A dt = 0 # ||dt|| <= sqrt(trust_radius**2 - ||dn||**2) # lb - dn <= dt <= ub - dn c_t = H . dot ( dn ) + c b_t = np . zeros_like ( b ) trust_radius_t = np . sqrt ( trust_radius ** 2 - np . linalg . norm ( dn ) ** 2 ) lb_t = trust_lb - dn ub_t = trust_ub - dn dt , cg_info = projected_cg ( H , c_t , Z , Y , b_t , trust_radius_t , lb_t , ub_t ) # Compute update (normal + tangential steps). d = dn + dt # Compute second order model: 1/2 d H d + c.T d + f. quadratic_model = 1 / 2 * ( H . dot ( d ) ) . dot ( d ) + c . T . dot ( d ) # Compute linearized constraint: l = A d + b. linearized_constr = A . dot ( d ) + b # Compute new penalty parameter according to formula (3.52), # reference [2]_, p.891. vpred = norm ( b ) - norm ( linearized_constr ) # Guarantee `vpred` always positive, # regardless of roundoff errors. vpred = max ( 1e-16 , vpred ) previous_penalty = penalty if quadratic_model > 0 : new_penalty = quadratic_model / ( ( 1 - PENALTY_FACTOR ) * vpred ) penalty = max ( penalty , new_penalty ) # Compute predicted reduction according to formula (3.52), # reference [2]_, p.891. predicted_reduction = - quadratic_model + penalty * vpred # Compute merit function at current point merit_function = f + penalty * norm ( b ) # Evaluate function and constraints at trial point x_next = x + S . dot ( d ) f_next , b_next = fun_and_constr ( x_next ) # Compute merit function at trial point merit_function_next = f_next + penalty * norm ( b_next ) # Compute actual reduction according to formula (3.54), # reference [2]_, p.892. actual_reduction = merit_function - merit_function_next # Compute reduction ratio reduction_ratio = actual_reduction / predicted_reduction # Second order correction (SOC), reference [2]_, p.892. if reduction_ratio < SUFFICIENT_REDUCTION_RATIO and norm ( dn ) <= SOC_THRESHOLD * norm ( dt ) : # Compute second order correction y = - Y . dot ( b_next ) # Make sure increment is inside box constraints _ , t , intersect = box_intersections ( d , y , trust_lb , trust_ub ) # Compute tentative point x_soc = x + S . dot ( d + t * y ) f_soc , b_soc = fun_and_constr ( x_soc ) # Recompute actual reduction merit_function_soc = f_soc + penalty * norm ( b_soc ) actual_reduction_soc = merit_function - merit_function_soc # Recompute reduction ratio reduction_ratio_soc = actual_reduction_soc / predicted_reduction if intersect and reduction_ratio_soc >= SUFFICIENT_REDUCTION_RATIO : x_next = x_soc f_next = f_soc b_next = b_soc reduction_ratio = reduction_ratio_soc # Readjust trust region step, formula (3.55), reference [2]_, p.892. if reduction_ratio >= LARGE_REDUCTION_RATIO : trust_radius = max ( TRUST_ENLARGEMENT_FACTOR_L * norm ( d ) , trust_radius ) elif reduction_ratio >= INTERMEDIARY_REDUCTION_RATIO : trust_radius = max ( TRUST_ENLARGEMENT_FACTOR_S * norm ( d ) , trust_radius ) # Reduce trust region step, according to reference [3]_, p.696. elif reduction_ratio < SUFFICIENT_REDUCTION_RATIO : trust_reduction = ( 1 - SUFFICIENT_REDUCTION_RATIO ) / ( 1 - reduction_ratio ) new_trust_radius = trust_reduction * norm ( d ) if new_trust_radius >= MAX_TRUST_REDUCTION * trust_radius : trust_radius *= MAX_TRUST_REDUCTION elif new_trust_radius >= MIN_TRUST_REDUCTION * trust_radius : trust_radius = new_trust_radius else : trust_radius *= MIN_TRUST_REDUCTION # Update iteration if reduction_ratio >= SUFFICIENT_REDUCTION_RATIO : x = x_next f , b = f_next , b_next c , A = grad_and_jac ( x ) S = scaling ( x ) # Get projections Z , LS , Y = projections ( A , factorization_method ) # Compute least-square lagrange multipliers v = - LS . dot ( c ) # Compute Hessian H = lagr_hess ( x , v ) # Set Flag last_iteration_failed = False # Otimality values optimality = norm ( c + A . T . dot ( v ) , np . inf ) constr_violation = norm ( b , np . inf ) if len ( b ) > 0 else 0 else : penalty = previous_penalty last_iteration_failed = True return x , state