Algorithms Homework 4

Course policies and AI category¶

Read the Artificial Intelligence Policy and Collaboration Policy and Honor Code before starting. Assignment-specific directions control. The categories are No AI, AI permitted after independent work, and AI required.

Unless a problem says otherwise, its category is AI permitted after independent work. Spend about 30 minutes on each top-level problem without AI, solution pages, or another person’s help, stopping early if complete. You may consult lecture notes, textbooks, and nonsolution pages of the course website; bias toward those course sources. Afterward, AI and genuine collaboration, including coding together, are permitted. Everyone must contribute intellectually, understand the work, and verify it.

At the end of each top-level problem, add a concise AI and independent-work report: approximately how long the independent attempt took, how far you got, where you became stuck, any AI or collaborative help used afterward, and how you verified it. If you used no AI, say so. Do not submit prompts or transcripts. Time estimates help the instructor improve the assignment and are not a speed test.

Assignment Overview: This assignment implements a basic interior point method for nonlinear programs with equality constraints and bounds.

Tips and Tricks¶

Background¶

Recall, the step for a primal-dual interior point method is defined by the following system of linear equations (Eq. 6.56 in Biegler, 2010):

\begin{bmatrix} W^k & \nabla c(x^k) & -I \\ \nabla c(x^k)^T & 0 & 0 \\ U^k & 0 & X^k \end{bmatrix} \begin{bmatrix} d_x^k \\ d_v^k \\ d_u^k \end{bmatrix} = -\begin{bmatrix} \nabla f(x^k) + \nabla c(x^k) v^k - u^k \\ c(x^k) \\ X^k u^k - \mu_l e \end{bmatrix}.

Alternatively, the step for the dual variables for the bounds can be recommed for this system (Eq. 6.57 in Biegler, 2010):

\begin{bmatrix} W^k + \Sigma^k & \nabla c(x^k) \\ \nabla c(x^k)^T & 0 \end{bmatrix} \begin{bmatrix} d_x^k \\ d_v^k \end{bmatrix} = -\begin{bmatrix} \nabla \varphi_\mu(x^k) + \nabla c(x^k) v^k \\ c(x^k) \end{bmatrix}.

And then computed using the solution to the linear system (Eq. 6.58 in Biegler, 2010):

d_u^k = \mu_l (X^k)^{-1} e - u^k - \Sigma^k d_x^k.

Finally, inertia correction and regularization can be applied to ensure reliable calculation of the Newton step (Eq. 6.59 in Biegler, 2010):

\begin{bmatrix} W^k + \Sigma^k + \delta_W I & \nabla c(x^k) \\ \nabla c(x^k)^T & -\delta_A I \end{bmatrix} \begin{bmatrix} d_x^k \\ d_v^k \end{bmatrix} = -\begin{bmatrix} \nabla \varphi_\mu(x^k) + \nabla c(x^k) v^k \\ c(x^k) \end{bmatrix}.

Problem Formulation¶

Consider the following nonlinear program:

\begin{align*} \min_{x} \quad & f(x) \\ \mathrm{s.t.} \quad & c(x) = 0 \\ & x_i \geq 0, \quad i \in \mathcal{I} \end{align*}

(1)

where $x \in \mathbb{R}^{n}$ , $f(x): \mathbb{R}^{n} \rightarrow \mathbb{R}$ , $c(x): \mathbb{R}^{n} \rightarrow \mathbb{R}^{m}$ and $|\mathcal{I} |$ = $r$ (i.e., there are $r$ variables with a lower bound and $r \leq n$ ). This is an extension of (6.48) in Biegler (2010).

This has the corresponding log-barrier approximation:

\begin{align*} \min_{x} \quad & \phi_{\mu_l}(x) := f(x) - \mu_l \sum_{i \in \mathcal{I}} \log(x_i) \\ \mathrm{s.t.} \quad & c(x) = 0 \end{align*}

(2)

which is an extension of (6.49) in Biegler (2010).

Let the $n \times r$ matrix $G$ encode which variables are bounded. If variable $i$ corresponds to the $j$ th bound, then $G_{i,j} = 1$ and otherwise $G_{i,j} = 0$ .

$G$ is assembled as follows:

Initialize $G$ as the zero matrix
Loop over $k=1$ to $k=|\mathcal{I}|$
1. Extract the index $i$ corresponding to the element $\mathcal{I}_k$
2. Set $G_{i,k} = 1$

Notice that $G$ is the gradient of $x \geq 0$ . $G^T G = I$ but the converse does not hold unless $r = n$ .

Reformulation Example¶

Start with:

\begin{align*} \min_{x} \quad & x_1^2 + x_2 \\ \mathrm{s.t.} \quad & x_1 + x_2 = 1 \\ & x_1 + 1 \geq 0 \end{align*}

(3)

Add slack variable $x_3$ and convert the inequality constraint to an equality constraint and bound:

\begin{align*} \min_{x} \quad & x_1^2 + x_2 \\ \mathrm{s.t.} \quad & x_1 + x_2 = 1 \\ & x_1 + 1 - x_3 = 0 \\ & x_3 \geq 0 \end{align*}

(4)

Now assemble $G$ :

G = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}

Primal Dual Optimality Conditions¶

Next we extend (6.51) in Biegler (2010):

\begin{gather*} \underbrace{\nabla f(x^k)}_{n \times 1} ~+~ \underbrace{\nabla c(x^k)}_{n \times m} \underbrace{v^k}_{m \times 1} ~-~ \underbrace{G}_{n \times r} \underbrace{u^k}_{r \times 1} = \underbrace{0}_{n \times 1} \\ \underbrace{\hat{X}^k}_{r \times r} \underbrace{u^k}_{r \times 1} = \mu \underbrace{e}_{r \times 1} \\ \underbrace{c(x^k)}_{m \times 1} = \underbrace{0}_{m \times 1} \end{gather*}

(5)

We now extend (6.56) in Biegler (2010):

\begin{bmatrix} \underbrace{W^k}_{n \times n} & \underbrace{\nabla c(x^k)}_{n \times m} & -\underbrace{G}_{n \times r} \\ \underbrace{\nabla_x c(x^k)^T}_{m \times n} & \underbrace{0}_{m \times m} & \underbrace{0}_{m \times r} \\ \underbrace{U^k}_{r \times r} \underbrace{G^T}_{r \times n} & \underbrace{0}_{r \times m} & \underbrace{\hat{X}^k}_{r \times r} \end{bmatrix} \begin{bmatrix} \underbrace{d_x^k}_{n \times 1} \\ \underbrace{d_v^k}_{m \times 1} \\ \underbrace{d_u^k}_{r \times 1} \end{bmatrix} = - \begin{bmatrix} \underbrace{\nabla f(x^k)}_{n \times 1} ~+~ \underbrace{\nabla c(x^k)}_{n \times m} \underbrace{v^k}_{m \times 1} ~-~ \underbrace{G}_{n \times r} \underbrace{u^k}_{r \times 1} \\ \underbrace{c(x^k)}_{m \times 1} \\ \underbrace{\hat{X}^k}_{r \times r} \underbrace{u^k}_{r \times 1} ~-~ \mu_l \underbrace{e}_{r \times 1} \end{bmatrix}

where $U^k = \mathrm{diag}\{u^k\}$ and $W^k = \nabla_{xx} L(x^k,v^k)$ . Notice that $W^k$ does NOT include a contribution from the barrier term:

\underbrace{W^k}_{n \times n} = \underbrace{\nabla^2 f(x^k)}_{n \times n} + \sum_{j=1}^{m} \left( \underbrace{\nabla^2 c_{j}(x^k)}_{n \times n} \underbrace{v^k_{j}}_{\mathrm{scalar}} \right)

We can verify that $W^k$ does not include $\hat{X}$ by showing that the KKT system above is a Newton step to solve the nonlinear system for the primal dual conditions.

The Newton step can be simplified, similar to (6.57) and (6.58) in Biegler (2010):

\begin{bmatrix} \underbrace{W^k}_{n \times n} + \underbrace{\Sigma^k}_{n \times n} & \underbrace{\nabla c(x^k)}_{n \times m} \\ \underbrace{\nabla_x c(x^k)^T}_{m \times n} & 0 \end{bmatrix} \begin{bmatrix} \underbrace{d_x^k}_{n \times 1} \\ \underbrace{d_v^k}_{m \times 1} \end{bmatrix} = - \begin{bmatrix} \underbrace{\nabla \phi_{\mu_l}}_{n \times 1} ~+~ \underbrace{\nabla c(x^k)}_{n \times m} \underbrace{v^k}_{m \times 1} \\ \underbrace{c(x^k)}_{m \times 1} \end{bmatrix}

and

\underbrace{d_u^k}_{r \times 1} = \mu_l \underbrace{(\hat{X}^k)^{-1}}_{r \times r} \underbrace{e}_{r \times 1} ~-~ \underbrace{u^k}_{r \times 1} ~-~ \underbrace{G^T}_{r \times n} \underbrace{\Sigma^k}_{n \times n} \underbrace{d_x^k}_{n \times 1}

where

\underbrace{\Sigma^k}_{n \times n} = \underbrace{G}_{n \times r} \underbrace{(\hat{X}^k)^{-1}}_{r \times r} \underbrace{U^k}_{r \times r} \underbrace{G^T}_{r \times n}

and

\underbrace{\nabla \phi_{\mu_l}}_{n \times 1} = \underbrace{\nabla f(x^k)}_{n \times 1} - \mu_l \underbrace{G}_{n \times r} \underbrace{(\hat{X}^k)^{-1}}_{r \times r} \underbrace{e}_{r \times 1}

Notice that the equation for $\underbrace{d_u^k}_{r \times 1}$ simplifies by substituting $G^T G = I$ :

\underbrace{d_u^k}_{r \times 1} = \mu_l \underbrace{(\hat{X}^k)^{-1}}_{r \times r} \underbrace{e}_{r \times 1} ~-~ \underbrace{u^k}_{r \times 1} ~-~ \underbrace{(\hat{X}^k)^{-1}}_{r \times r} \underbrace{U^k}_{r \times r} \underbrace{G^T}_{r \times n} \underbrace{d_x^k}_{n \times 1}

Finally, inertia correction can be applied to simplified KKT step similar to (6.59) in Biegler (2010).

Globalization¶

Everything above computes a direction $d_x^k$ . Newton’s method is only locally convergent: from a starting point far from the solution, the full step $x^{k+1} = x^k + d_x^k$ can leave the interior, increase the objective, or diverge outright. A globalization strategy decides how far along $d_x^k$ to actually move.

We will implement and compare three strategies. All three share the same step calculation; they differ only in how the step length $\alpha^k$ is chosen.

Strategy 0: none. Always take the full step, $\alpha^k = 1$ . This is Part 1 above.

Strategy 1: $\ell_1$ merit function. Backtrack until a single scalar trade-off between the objective and the constraint violation decreases.

Strategy 2: filter. Backtrack until the pair (constraint violation, objective) is not dominated by anything seen before.

Strategies 1 and 2 both start their backtracking from the largest step that keeps the iterate strictly interior --- the step-to-the-boundary rule below.

The implementation also carries a fourth setting, "boundary", which applies the step-to-the-boundary rule and then accepts that step with no test at all. It is not a globalization strategy; it exists so we can measure how much of the improvement comes from the merit function or the filter and how much comes from merely staying interior.

Step-to-the-boundary rule¶

The log barrier $-\mu_l \sum_{i \in \mathcal{I}} \log(x_i)$ is only defined for $x_i > 0$ . Before any acceptance test we cap the step so the iterate stays strictly inside the bounds (Eq. 6.61a and 6.62 in Biegler, 2010):

\alpha_{\max}^k := \max \left\{ \alpha \in (0,1] ~:~ x^k + \alpha \, d_x^k \geq (1 - \tau_l) \, x^k \right\}, \qquad \tau_l = \max\{\tau_{\min}, 1 - \mu_l\},

with $\tau_{\min} \approx 0.99$ . Because $G^T x^k > 0$ componentwise, this is a scalar minimum over the bounded components with a negative step:

\alpha_{\max}^k = \min \left( 1, ~ \min_{i \in \mathcal{I} ~:~ (d_x^k)_i < 0} \frac{-\tau_l \, x_i^k}{(d_x^k)_i} \right).

The bound multipliers $u^k$ get their own step length $\alpha_u^k$ from the same rule applied to $u^k$ and $d_u^k$ (Eq. 6.61b, used in Eq. 6.60c in Biegler, 2010). Note $\tau_l \to 1$ as $\mu_l \to 0$ , so late iterations are allowed to come very close to the boundary.

Strategy 1: the $\ell_1$ merit function¶

Inside a barrier subproblem the “objective” is the barrier objective and the “infeasibility” is measured in the $\ell_1$ norm:

\varphi_{\mu_l}(x) = f(x) - \mu_l \sum_{i \in \mathcal{I}} \log(x_i), \qquad \theta(x) = \| c(x) \|_1 .

The $\ell_1$ merit function combines them with a penalty parameter $\rho > 0$ (Eq. 5.59, pg. 111 in Biegler, 2010):

\phi_1(x; \rho) = \varphi_{\mu_l}(x) + \rho \, \| c(x) \|_1 .

$\phi_1$ is nonsmooth, so the Armijo test uses the directional derivative rather than a gradient (Eq. 5.67, pg. 115 in Biegler, 2010):

D_{d_x} \phi_1(x^k; \rho) = \nabla \varphi_{\mu_l}(x^k)^T d_x^k - \rho \, \| c(x^k) \|_1 .

Only the second term is reliably negative. $\rho$ is not a tuning knob --- it is what makes $d_x^k$ a descent direction at all. It must exceed $\| v^* \|_\infty$ , which we do not know, so we use the multiplier estimate the Newton step just produced (Eq. 5.74, pg. 117 in Biegler, 2010):

\rho^k = \max \left\{ \rho^{k-1}, ~ \| v^k + d_v^k \|_\infty + \varrho \right\}, \qquad \varrho > 0 .

The line search is then the ordinary backtracking loop of Algorithm 3.2 with $\phi_1$ substituted for $f$ (Algorithm 5.3, pg. 117--118 in Biegler, 2010). Accept the first $\alpha \in \{\alpha_{\max}^k, \, \xi \alpha_{\max}^k, \, \xi^2 \alpha_{\max}^k, \dots\}$ satisfying

\phi_1(x^k + \alpha d_x^k; \rho^k) \leq \phi_1(x^k; \rho^k) + \eta \, \alpha \, D_{d_x}\phi_1(x^k; \rho^k), \qquad \eta \in (0, \tfrac{1}{2}) .

Strategy 2: the filter¶

The filter refuses to pick a trade-off weight at all. Read the barrier subproblem as a biobjective problem --- minimize $\varphi_{\mu_l}(x)$ and minimize $\theta(x)$ --- and keep a list $\mathcal{F}$ of the $(\theta, \varphi_{\mu_l})$ pairs seen so far.

A pair $(\theta_k, \varphi_k)$ dominates $(\theta_l, \varphi_l)$ if $\theta_k \leq \theta_l$ and $\varphi_k \leq \varphi_l$ . A trial point is acceptable to the filter if it is dominated by no entry of $\mathcal{F}$ .

Pure non-dominance is not quite enough --- a trial point could creep along the boundary of the forbidden region forever --- so acceptance carries a margin proportional to the current infeasibility (Eq. 5.60, pg. 113 in Biegler, 2010). A trial point $\hat{x}$ is accepted if

\theta(\hat{x}) \leq (1 - \gamma_\theta) \, \theta(x^k) \qquad \text{or} \qquad \varphi_{\mu_l}(\hat{x}) \leq \varphi_{\mu_l}(x^k) - \gamma_\varphi \, \theta(x^k) ,

for small $\gamma_\theta, \gamma_\varphi \in (0,1)$ .

Why a switching condition is needed. Once an iterate is feasible, $\theta(x^k) = 0$ and the test above degenerates: the second branch demands $\varphi_{\mu_l}(\hat{x}) \leq \varphi_{\mu_l}(x^k) - 0$ , i.e. essentially no decrease at all. The filter would happily converge to a feasible but non-optimal point. So near feasibility the filter has to stop being a filter and start being an Armijo line search on $\varphi_{\mu_l}$ alone.

That is what the switching condition does. Let $m_k(\alpha) := \alpha \, \nabla \varphi_{\mu_l}(x^k)^T d_x^k$ be the predicted decrease. When $m_k(\alpha) < 0$ and

\left[ -m_k(\alpha) \right]^{s_\varphi} \left[ \alpha \right]^{1 - s_\varphi} > \delta \left[ \theta(x^k) \right]^{s_\theta}, \qquad s_\theta > 1, ~ s_\varphi \geq 1 ,

the algorithm abandons the filter test and demands the ordinary Armijo condition on the barrier objective alone,

\varphi_{\mu_l}(\hat{x}) \leq \varphi_{\mu_l}(x^k) + \eta \, m_k(\alpha) ,

and does not augment the filter --- $\varphi_{\mu_l}$ strictly decreased, and that alone prevents cycling. See Section 5.6.2, pg. 119--120 and Algorithm 5.4, pg. 120--121 in Biegler (2010).

Part 1: Basic Interior Point Method for Inequality and Equality Constraint NLPs¶

Implement a basic interior point method for inequality and equality constrained nonlinear programs. See pg. 154 – 155 in Biegler (2010). You may skip the line search, i.e., always take a full step.

In Part 2 you will add the line search back and compare three ways of doing it.

Pseudocode¶

Write detailed pseudocode on paper or a whiteboard. Scan/take a photo and turn in.

Python Implementation¶

Implement in Python. Hints: Reuse code from Algorithm 5.2 example.

### Load Python libraries

import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper

    # helper.easy_install() # We do NOT need Pyomo for this assignment
else:
    sys.path.insert(0, "../../../optimization/notebooks/")
    import helper
helper.set_plotting_style()  # But we do want the nice plots

import numpy as np
from scipy import linalg


### Define helper functions
## Check is element of array is NaN
def check_nan(A):
    return np.sum(np.isnan(A))


## Calculate gradient with central finite difference
def my_grad_approx(x, f, eps1, verbose=False):
    """
    Calculate gradient of function f using central difference formula

    Inputs:
        x - point for which to evaluate gradient
        f - function to consider
        eps1 - perturbation size

    Outputs:
        grad - gradient (vector)
    """

    n = len(x)
    grad = np.zeros(n)

    if verbose:
        print("***** my_grad_approx at x = ", x, "*****")

    for i in range(0, n):

        # Create vector of zeros except eps in position i
        e = np.zeros(n)
        e[i] = eps1

        # Finite difference formula
        my_f_plus = f(x + e)
        my_f_minus = f(x - e)

        # Diagnostics
        if verbose:
            print("e[", i, "] = ", e)
            print("f(x + e[", i, "]) = ", my_f_plus)
            print("f(x - e[", i, "]) = ", my_f_minus)

        grad[i] = (my_f_plus - my_f_minus) / (2 * eps1)

    if verbose:
        print("***** Done. ***** \n")

    return grad


## Calculate gradient with central finite difference
def my_jac_approx(x, h, eps1, verbose=False):
    """
    Calculate Jacobian of function h(x) using central difference formula

    Inputs:
        x - point for which to evaluate gradient
        h - vector-valued function to consider. h(x): R^n --> R^m
        eps1 - perturbation size

    Outputs:
        A - Jacobian (n x m matrix)
    """

    # Check h(x) at x
    h_x0 = h(x)

    # Extract dimensions
    n = len(x)
    m = len(h_x0)

    # Initialize Jacobian matrix
    A = np.zeros((n, m))

    # Calculate Jacobian by row
    for i in range(0, n):

        # Create vector of zeros except eps in position i
        e = np.zeros(n)
        e[i] = eps1

        # Finite difference formula
        my_h_plus = h(x + e)
        my_h_minus = h(x - e)

        # Diagnostics
        if verbose:
            print("e[", i, "] = ", e)
            print("h(x + e[", i, "]) = ", my_h_plus)
            print("h(x - e[", i, "]) = ", my_h_minus)

        A[i, :] = (my_h_plus - my_h_minus) / (2 * eps1)

    if verbose:
        print("***** Done. ***** \n")

    return A


## Calculate Hessian using central finite difference
def my_hes_approx(x, grad, eps2):
    """
    Calculate gradient of function my_f using central difference formula and my_grad

    Inputs:
        x - point for which to evaluate gradient
        grad - function to calculate the gradient
        eps2 - perturbation size (for Hessian NOT gradient approximation)

    Outputs:
        H - Hessian (matrix)
    """

    n = len(x)
    H = np.zeros([n, n])

    for i in range(0, n):
        # Create vector of zeros except eps in position i
        e = np.zeros(n)
        e[i] = eps2

        # Evaluate gradient twice
        grad_plus = grad(x + e)
        grad_minus = grad(x - e)

        # Notice we are building the Hessian by column (or row)
        H[:, i] = (grad_plus - grad_minus) / (2 * eps2)

    return H


## Linear algebra calculation
def xxT(u):
    """
    Calculates u*u.T to circumvent limitation with SciPy

    Arguments:
    u - numpy 1D array

    Returns:
    u*u.T

    Assume u is a nx1 vector.
    Recall: NumPy does not distinguish between row or column vectors

    u.dot(u) returns a scalar. This functon returns an nxn matrix.
    """

    n = len(u)
    A = np.zeros([n, n])
    for i in range(0, n):
        for j in range(0, n):
            A[i, j] = u[i] * u[j]

    return A


## Analyze Hessian
def analyze_hes(B):
    print(B, "\n")

    l = linalg.eigvals(B)
    print("Eigenvalues: ", l, "\n")

## Assemble KKT matrix (equality constrained only)
def assemble_check_KKT(W, Sk, A, deltaA, deltaW, verbose):

    # Add your solution here

    return KKT, inertia_correct, pos_ev, neg_ev, zero_ev

Adding globalization to the implementation¶

We now extend barrier_subproblem with a globalization keyword argument taking the values "none", "boundary", "merit", or "filter". The default is "none", so all of the Part 1 results above are unchanged.

Both barrier_subproblem and interior_point also accept an optional info dictionary. If you pass one in, it is filled with per-iteration diagnostics --- step lengths, $\theta$ , $\varphi_{\mu_l}$ , $\rho$ , the filter, the number of trial points --- and a termination status. We use it below to build the comparison.

Start with the four short helper functions.

### Globalization helper functions


def barrier_objective(x, calc_f, var_bounds, mu):
    """
    Evaluate the log-barrier objective phi_mu(x) = f(x) - mu * sum(log(x_i))

    Inputs:
        x - point at which to evaluate (vector)
        calc_f - function to calculate objective (returns scalar)
        var_bounds - list of indices for variables with a lower bound of zero
        mu - barrier penalty

    Outputs:
        phi - barrier objective. Returns np.inf if x is outside the domain
              of the log, which lets a line search reject the trial point
              instead of crashing.
    """

    # Add your solution here

    return phi


def constraint_violation(x, calc_c):
    """
    Evaluate the infeasibility measure theta(x) = ||c(x)||_1

    Inputs:
        x - point at which to evaluate (vector)
        calc_c - function to calculate constraints (returns vector)

    Outputs:
        theta - 1-norm of the constraint residual (scalar)
    """

    # Add your solution here

    return theta


def fraction_to_boundary(z, dz, tau):
    """
    Step-to-the-boundary rule, Eq. 6.61a and 6.62 in Biegler (2010)

    Largest alpha in (0, 1] such that z + alpha*dz >= (1 - tau)*z, elementwise.

    Inputs:
        z - current value of a strictly positive vector (x on the bounded
            components, or the bound multipliers u)
        dz - Newton step for z
        tau - fraction-to-boundary parameter, max(tau_min, 1 - mu)

    Outputs:
        alpha_max - largest permissible step length (scalar in (0, 1])
    """

    # Add your solution here

    return alpha_max


def filter_dominated(theta, phi, filter_set):
    """
    Is the pair (theta, phi) dominated by any entry in the filter?

    A filter entry (theta_j, phi_j) dominates (theta, phi) if it is no worse
    in BOTH coordinates. A dominated trial point is unacceptable.

    Inputs:
        theta - constraint violation of the trial point (scalar)
        phi - barrier objective of the trial point (scalar)
        filter_set - list of (theta_j, phi_j) tuples

    Outputs:
        True if (theta, phi) is dominated, False otherwise
    """

    # Add your solution here

    return False

def barrier_subproblem(
    x0,
    v0,
    u0,
    calc_f,
    calc_c,
    var_bounds,
    mu,
    eps,
    max_iter=100,
    verbose=False,
    globalization="none",
    info=None,
):
    """
    Basic Full Space Newton Method for Solving Barrier Subproblem Constrained NLP

    Input:
        x0 - starting point (vector)
        calc_f - function to calculate objective (returns scalar)
        calc_c - function to calculate constraints (returns vector)
        var_bounds - list of indices for variables with lower bound
        mu - barrier penalty
        eps - tolerance for termination
        globalization - step acceptance strategy. One of:
            "none"     - always take the full step, alpha = 1 (Part 1)
            "boundary" - step-to-the-boundary rule only, no acceptance test
            "merit"    - backtracking line search on the l1 merit function
            "filter"   - backtracking filter line search
        info - optional dictionary. If provided, it is filled with per-iteration
               diagnostics and a termination status.

    Histories (stored for debugging):
        x - history of steps (primal variables)
        v - history of steps (duals for constraints)
        u - history of steps (duals for bounds)
        f - history of objective evaluations
        c - history of constraint evaluations
        df - history of objective gradients
        dL - history of Lagrange function gradients
        A - history of constraint Jacobians
        W - history of Lagrange Hessians
        S - history of sigma matrix

    Outputs:
        x - final value for primal variable
        v - final value for constraint duals
        u - final value for bound duals

    Notes:
        1. For simplicity, central finite difference is used
           for all gradient calculations.
    """

    ### Specifics for Algorithm 5.2
    # Tuning parameters
    delta_bar_W_min = 1e-20
    delta_bar_W_0 = 1e-4
    delta_bar_W_max = 1e40
    delta_bar_A = 1e-8
    kappa_u = 8
    kappa_l = 1 / 3

    ### Specifics for globalization
    # Step-to-the-boundary rule, Eq. 6.62 in Biegler (2010)
    tau_min = 0.99
    # Armijo constant, shared by the merit function and the filter
    eta_armijo = 1e-4
    # Backtracking factor and smallest step length we will try
    xi_backtrack = 0.5
    alpha_min = 1e-10
    # Safety margin in the penalty parameter update, Eq. 5.74 in Biegler (2010)
    rho_margin = 1.0
    # Filter margins, Eq. 5.60 in Biegler (2010)
    gamma_theta = 1e-5
    gamma_phi = 1e-5
    # Switching condition, Section 5.6.2 in Biegler (2010)
    delta_switch = 1.0
    s_phi = 2.3
    s_theta = 1.1

    if globalization not in ("none", "boundary", "merit", "filter"):
        raise ValueError("Unknown globalization strategy: " + str(globalization))

    # Declare iteration histories as empty lists
    x = []
    v = []
    u = []
    f = []
    L = []
    c = []
    df = []
    dL = []
    A = []
    W = []

    # Add your solution here

def interior_point(
    x0, calc_f, calc_c, var_bounds, max_iter=20, globalization="none", info=None
):
    """
    Basic interior point method: a sequence of barrier subproblems

    Inputs:
        x0 - starting point (vector)
        calc_f - function to calculate objective (returns scalar)
        calc_c - function to calculate constraints (returns vector)
        var_bounds - list of indices for variables with lower bound
        max_iter - maximum number of barrier subproblems
        globalization - "none", "boundary", "merit", or "filter"
        info - optional dictionary filled with diagnostics

    Outputs:
        x, v, u, mu, E - histories over the barrier subproblems
    """

    # Add your solution here

    if info is not None:
        info["status"] = status
        info["iterations"] = all_iterations
        info["filter"] = filter_set
        info["filters"] = filters
        info["rho"] = rho
        info["newton_iterations"] = len(all_iterations)
        info["trial_points"] = sum(e["trials"] for e in all_iterations)
        info["barrier_subproblems"] = k

    return x, v, u, mu, E

Test Problems¶

Problem 1: Convex¶

f = lambda x: x[0] + 2 * x[1]
c = lambda x: (x[0] + x[1] - 1) * np.ones(1)

# Indices of variables with lower bound of zero
vb = [1]

x0 = np.ones(2)
u0 = np.ones(1)
v0 = np.ones(1)

x_, v_, u_, E_ = barrier_subproblem(x0, v0, u0, f, c, vb, 1e-1, 1e-10, verbose=False)

Problem 2: Convex¶

f = lambda x: x[0] + 2 * x[1]
c = lambda x: (x[0] + x[1] - 1) * np.ones(1)

# Indices of variables with lower bound of zero
vb = [1]

x0 = np.ones(2)
# u0 = np.ones(1)
# v0 = np.ones(1)

# x_, v_, u_, E_ = barrier_subproblem(x0,v0,u0,f,c,vb,1E-1,1E-10,verbose=False)

x, v, u, mu, E = interior_point(x0, f, c, vb)

import matplotlib.pyplot as plt


def plot_results(x, u, v, mu, E):
    """Plot results of interior point method

    Inputs:
        x - primal variables
        u - dual variables for lower bounds
        v - dual variables for equality constraints
        mu - barrier penalty
        E - error

    Outputs:
        None

    Other:
        Plots are generated.
    """

    # number of iterations
    N = len(x)

    # iteration
    iters = range(0, N)

    plt.figure()
    for i in range(0, len(x[0])):
        plt.plot(iters, [x[j][i] for j in range(0, N)], label="$x_{" + str(i) + "}$")
    plt.xlabel("Iteration")
    plt.ylabel("Primal Variables")
    plt.legend()
    plt.show()

    plt.figure()
    for i in range(0, len(u[0])):
        plt.plot(iters, [u[j][i] for j in range(0, N)], label="$u_{" + str(i) + "}$")
    for i in range(0, len(v[0])):
        plt.plot(iters, [v[j][i] for j in range(0, N)], label="$v_{" + str(i) + "}$")
    plt.xlabel("Iteration")
    plt.ylabel("Dual Variables")
    plt.legend()
    plt.show()

    plt.figure()
    plt.semilogy(iters, mu)
    plt.xlabel("Iteration")
    plt.ylabel("Barrier Penalty")
    plt.show()

    plt.figure()
    plt.semilogy(range(1, N), E)
    plt.xlabel("Iteration")
    plt.ylabel("Error")
    plt.show()


plot_results(x, u, v, mu, E)

Problem 3: Nonconvex¶

f = lambda x: x[0] + 2 * x[1] + x[2] ** 2


def c(x):
    rhs = np.zeros(2)
    rhs[0] = x[0] + x[1] - 1
    rhs[1] = x[2] - x[1] - 3
    return rhs


# Indices of variables with lower bound of zero
vb = [1, 2]

x0 = np.ones(3)
u0 = np.ones(2)
v0 = np.ones(2)

## Test barrier subproblem
# xtest, vtest, utest, Etest = barrier_subproblem(x0,v0,u0,f,c,vb,1E-1,1E-10,max_iter=10,verbose=True)

##
x, v, u, mu, E = interior_point(x0, f, c, vb)

plot_results(x, u, v, mu, E)

Comparing Globalization Strategies¶

Problems 1--3 all converge with no globalization at all. That is not an accident --- they were chosen in Part 1 so the bare Newton step works. The interesting question is what happens when it does not.

The three problems below are chosen so the strategies actually separate. Each one is solved three times --- "none", "merit", "filter" --- from the same starting point, and every answer is cross-checked against an analytic solution and against Ipopt.

### Setup for the comparison

import contextlib
import io
import warnings

# Okabe-Ito colorblind-safe palette (see figures/README.md in the public repo)
OKABE_ITO = {
    "blue": "#0072B2",
    "orange": "#E69F00",
    "vermillion": "#D55E00",
    "bluishgreen": "#009E73",
    "skyblue": "#56B4E9",
    "reddishpurple": "#CC79A7",
    "yellow": "#F0E442",
    "black": "#000000",
}

# Line style per strategy, so every figure below is read the same way
STRATEGY_STYLE = {
    "none": (OKABE_ITO["vermillion"], "-"),
    "boundary": (OKABE_ITO["orange"], "-."),
    "merit": (OKABE_ITO["blue"], "--"),
    "filter": (OKABE_ITO["bluishgreen"], ":"),
}

def compare_globalization(
    x0,
    calc_f,
    calc_c,
    var_bounds,
    strategies=("none", "merit", "filter"),
    max_iter=20,
    reference=None,
):
    """
    Solve the same NLP once per globalization strategy and tabulate the result.

    Inputs:
        x0 - starting point (vector), shared by every strategy
        calc_f - function to calculate objective (returns scalar)
        calc_c - function to calculate constraints (returns vector)
        var_bounds - list of indices for variables with a lower bound of zero
        strategies - globalization strategies to compare
        max_iter - maximum number of barrier subproblems
        reference - optional known solution x*, used to report ||x - x*||

    Outputs:
        results - dict mapping strategy name to (x, v, u, mu, E, info)
    """

    results = {}

    print(
        "Strategy   \tStatus              \tNewton \tTrial \tf(x)        \t||c(x)||_1  \t||x - x*||"
    )
    print("-" * 108)

    for g in strategies:
        info = {}

        # Suppress the per-iteration log from the inner solver. A strategy that is
        # failing also emits LinAlgWarning from linalg.solve; that is a symptom, not
        # an error, so it is silenced here to keep the table legible.
        with contextlib.redirect_stdout(io.StringIO()), warnings.catch_warnings():
            warnings.simplefilter("ignore", linalg.LinAlgWarning)
            x, v, u, mu, E = interior_point(
                x0.copy(),
                calc_f,
                calc_c,
                var_bounds,
                max_iter=max_iter,
                globalization=g,
                info=info,
            )

        results[g] = (x, v, u, mu, E, info)

        xf = np.asarray(x[-1], dtype=float)

        if np.all(np.isfinite(xf)):
            f_str = "{0: 1.5e}".format(calc_f(xf))
            c_str = "{0: 1.3e}".format(constraint_violation(xf, calc_c))
            if reference is not None:
                d_str = "{0: 1.3e}".format(
                    linalg.norm(xf - np.asarray(reference), np.inf)
                )
            else:
                d_str = "    -----"
        else:
            f_str, c_str, d_str = "    -----", "    -----", "    -----"

        print(
            "{0:11s}\t{1:20s}\t{2:6d} \t{3:5d} \t{4:s} \t{5:s} \t{6:s}".format(
                g,
                info["status"],
                info["newton_iterations"],
                info["trial_points"],
                f_str,
                c_str,
                d_str,
            )
        )

    return results

def plot_filter_pairs(info, title="", mu=None):
    """Plot the (theta, phi) pairs one barrier subproblem visited.

    The barrier objective phi_mu(x) = f(x) - mu*sum(log(x_i)) depends on mu, so
    pairs from different barrier subproblems are NOT comparable and must not
    share an axis. This function therefore draws a single subproblem: by default
    the one that needed the most trial points, i.e. the one where the line search
    actually had to work.

    Inputs:
        info - diagnostics dictionary returned by interior_point(..., info=...)
        title - optional axes title
        mu - barrier parameter selecting which subproblem to draw. Default:
             the subproblem with the most trial points.

    Outputs:
        None. A figure is generated.
    """

    if mu is None:
        totals = {}
        for e in info["iterations"]:
            totals[e["mu"]] = totals.get(e["mu"], 0) + e["trials"]
        mu = max(totals, key=totals.get)

    steps = [e for e in info["iterations"] if e["mu"] == mu and np.isfinite(e["phi"])]

    accepted = np.array([[e["theta"], e["phi"]] for e in steps])
    rejected = np.array(
        [[t, p] for e in steps for (t, p) in e["rejected"] if np.isfinite(p)]
    )
    entries = np.array(info.get("filters", {}).get(mu, []))

    fig, ax = plt.subplots(figsize=(6, 4.5))

    if len(rejected) > 0:
        ax.plot(
            rejected[:, 0],
            rejected[:, 1],
            linestyle="none",
            marker="x",
            markersize=9,
            markeredgewidth=2,
            color=OKABE_ITO["vermillion"],
            label="rejected trial points",
        )

    if len(accepted) > 0:
        ax.plot(
            accepted[:, 0],
            accepted[:, 1],
            linestyle="-",
            marker="o",
            markersize=7,
            color=OKABE_ITO["blue"],
            label="accepted iterates",
        )

    if len(entries) > 0:
        ax.plot(
            entries[:, 0],
            entries[:, 1],
            linestyle="none",
            marker="s",
            markersize=11,
            markerfacecolor="none",
            markeredgewidth=2,
            color=OKABE_ITO["orange"],
            label="filter entries",
        )

    ax.set_xscale("symlog", linthresh=1e-10)
    ax.set_yscale("symlog", linthresh=1.0)
    ax.set_xlabel(r"$\theta(x) = \|c(x)\|_1$")
    ax.set_ylabel(r"$\varphi_{\mu_l}(x)$, at $\mu_l = $ " + "{0:g}".format(mu))
    if title:
        ax.set_title(title + r" ($\mu_l = $ " + "{0:g})".format(mu))
    ax.legend(loc="best", fontsize=11)
    ax.tick_params(top=True, right=True, direction="in")
    plt.tight_layout()
    plt.show()

Problem 4: Where the full Newton step diverges¶

\begin{align*} \min_{x} \quad & \sqrt{1 + x_1^2} + \tfrac{1}{2} x_1 + x_2 \\ \mathrm{s.t.} \quad & x_2 - x_3 - 1 = 0 \\ & x_2 \geq 0, \quad x_3 \geq 0 \end{align*}

(6)

from $x^0 = (2, 2, 1)$ .

The objective separates. In $x_1$ it is the classic Newton counterexample: for $g(t) = \sqrt{1+t^2}$ the Newton step is $-g'/g'' = -t(1+t^2)$ , which overshoots and diverges whenever $|t| > 1$ . Nothing about the bounds or the constraint is difficult here; the trouble is entirely in an unbounded variable, so the step-to-the-boundary rule cannot help.

The analytic solution is available. With $x_2 = x_3 + 1$ and $x_3 \geq 0$ , the objective is minimized at $x_3 = 0$ , and stationarity in $x_1$ gives $x_1/\sqrt{1+x_1^2} = -\tfrac{1}{2}$ , so

x^* = \left( -\tfrac{1}{\sqrt{3}}, ~ 1, ~ 0 \right) \approx (-0.5773503, ~ 1, ~ 0), \qquad f(x^*) = \frac{\sqrt{3}}{2} + 1 \approx 1.8660254 .

### Problem 4: sqrt(1 + x1^2) + x1/2 + x2

f4 = lambda x: np.sqrt(1 + x[0] ** 2) + 0.5 * x[0] + x[1]


def c4(x):
    rhs = np.zeros(1)
    rhs[0] = x[1] - x[2] - 1
    return rhs


# Indices of variables with lower bound of zero
vb4 = [1, 2]

x0_4 = np.array([2.0, 2.0, 1.0])

# Analytic solution, derived in the markdown cell above
x_star_4 = np.array([-1 / np.sqrt(3), 1.0, 0.0])
f_star_4 = np.sqrt(3) / 2 + 1

print("Analytic solution: x* =", x_star_4, " f(x*) =", f_star_4)
print()

results4 = compare_globalization(
    x0_4,
    f4,
    c4,
    vb4,
    strategies=("none", "boundary", "merit", "filter"),
    reference=x_star_4,
)

### Problem 4: what the un-globalized method actually does

# The full Newton step is unstable in x1. Look at the iterates it produced.
x_none = results4["none"][0]
info_none = results4["none"][5]

print("globalization = 'none', status:", info_none["status"])
print("final iterate:", np.asarray(x_none[-1], dtype=float))
print(
    "objective at final iterate: {0:.4f}  (optimal value is {1:.4f})".format(
        f4(np.asarray(x_none[-1], dtype=float)), f_star_4
    )
)
print()
print("Iterate at the end of each barrier subproblem:")
for j in range(0, len(x_none)):
    print("  subproblem {0}: x = {1}".format(j, np.asarray(x_none[j], dtype=float)))
print()
print(
    "x1 is unbounded below and does not appear in c(x). The step-to-the-boundary\n"
    "rule cannot help there, which is why the 'boundary' row above also fails."
)

### Problem 4: convergence of the four settings

fig, ax = plt.subplots(figsize=(6.4, 4.5))

for g in ("none", "boundary", "merit", "filter"):
    info = results4[g][5]
    err = [e["E"] for e in info["iterations"]]
    err = [max(e, 1e-16) for e in err]
    color, dash = STRATEGY_STYLE[g]
    ax.semilogy(
        range(1, len(err) + 1),
        err,
        dash,
        color=color,
        label="{0} ({1})".format(g, info["status"]),
    )

ax.set_xlabel("Newton iteration")
ax.set_ylabel(r"KKT error $E$")
ax.set_title("Problem 4: KKT error per Newton iteration")
ax.tick_params(top=True, right=True, direction="in")
ax.legend(loc="best", fontsize=11)
plt.tight_layout()
plt.show()

### Problem 4: the filter line search, drawn

# Problem 5's filter run never rejects a single trial point, so there is nothing
# to draw there. Problem 4 is where the line search has to work: look at the one
# barrier subproblem that needed the most trial points.
#
# Read the x-axis carefully. c(x) here is LINEAR and x0 is already feasible, so
# theta is at roundoff (order 1e-10) for the whole run and the horizontal spread
# is floating-point noise, not structure. What the figure actually shows is the
# backtracking Armijo branch: nine trial points marching down in phi_mu until one
# finally drops below the current iterate. Neither of these test problems ever
# builds a filter with more than one entry -- the switching condition sends every
# step down the Armijo branch. That is the behavior the theory predicts near
# feasibility (Section 5.6.2 in Biegler, 2010), not a bug in our implementation.

plot_filter_pairs(results4["filter"][5], title=r"Problem 4: filter line search")

Problem 5: Where the merit function and the filter disagree¶

\begin{align*} \min_{x} \quad & 2 x_1 + x_2 \\ \mathrm{s.t.} \quad & x_1^2 + x_2^2 - 1 = 0 \\ & x_2 \geq 0 \end{align*}

(7)

from the infeasible starting point $x^0 = (0.9, 0.9)$ , where $\theta(x^0) = 0.62$ .

The constraint is a circle, so the linearization is poor when far from it. The unconstrained minimizer of $2x_1 + x_2$ on the circle is $-(2,1)/\sqrt{5}$ , which has $x_2 < 0$ ; the bound is therefore active and

x^* = (-1, ~ 0), \qquad f(x^*) = -2 .

This is the geometry of the acceptance regions, made concrete. At the very first iteration ( $\mu_0 = 10$ , $x^0 = (0.9, 0.9)$ ) the full step goes to $x^0 + d_x^0 = (-1.679, 3.135)$ and trades a large increase in infeasibility for a large decrease in the barrier objective:

	$\theta$	$\varphi_{\mu_0}$
$x^0$	0.620	3.754
$x^0 + d_x^0$	11.645	-11.648

The two strategies read that trade differently.

The merit function collapses it to one number. Here $\rho^0 = 2.755$ , so the step changes $\phi_1 = \varphi_{\mu_0} + \rho\,\theta$ by $(-15.40) + 2.755 \times (+11.02) = +14.97$ : the penalized infeasibility swamps the objective decrease, the Armijo test fails, and the search backtracks to $\alpha^0 = 0.5$ .

The filter never forms that sum. The predicted decrease is $m_0 = \nabla \varphi_{\mu_0}^T d_x^0 = -27.75$ , which is large enough relative to $\theta(x^0)^{s_\theta}$ that the switching condition (5.85) fires. The step is therefore an $f$ -type step: the filter asks only for Armijo decrease in $\varphi_{\mu_0}$ , gets it easily, accepts $\alpha^0 = 1$ , and does not augment the filter.

### Problem 5: 2*x1 + x2 on the unit circle

f5 = lambda x: 2 * x[0] + x[1]


def c5(x):
    rhs = np.zeros(1)
    rhs[0] = x[0] ** 2 + x[1] ** 2 - 1
    return rhs


# Indices of variables with lower bound of zero
vb5 = [1]

x0_5 = np.array([0.9, 0.9])

# Analytic solution
x_star_5 = np.array([-1.0, 0.0])
f_star_5 = -2.0

print("Analytic solution: x* =", x_star_5, " f(x*) =", f_star_5)
print("theta(x0) = {0:.4f}".format(constraint_violation(x0_5, c5)))
print()

results5 = compare_globalization(x0_5, f5, c5, vb5, reference=x_star_5)

print()
print("Step lengths at each Newton iteration:")
for g in ("none", "merit", "filter"):
    alphas = [round(float(e["alpha"]), 4) for e in results5[g][5]["iterations"]]
    print("  {0:7s} {1}".format(g, alphas))

print()
print("Penalty parameter rho (merit only): {0:.4f}".format(results5["merit"][5]["rho"]))

Problem 6: Where all three fail¶

Globalization is not a cure-all. This is the Wachter--Biegler counterexample, the problem that motivated the feasibility restoration phase. It is Example 6.8 (Eq. 6.63), pg. 155--156 in Biegler (2010), and we use the same starting point the book does:

\begin{align*} \min_{x} \quad & x_1 \\ \mathrm{s.t.} \quad & x_1^2 - x_2 - 1 = 0 \\ & x_1 - x_3 - \tfrac{1}{2} = 0 \\ & x_2 \geq 0, \quad x_3 \geq 0 \end{align*}

(8)

from $x^0 = (-2, 3, 1)$ .

The feasible set requires $x_1 \geq \tfrac{1}{2}$ (from the second constraint and $x_3 \geq 0$ ) and $x_1^2 \geq 1$ (from the first and $x_2 \geq 0$ ), so $x_1 \geq 1$ and

x^* = (1, ~ 0, ~ \tfrac{1}{2}), \qquad f(x^*) = 1 .

The solution exists and is perfectly well behaved. But from a starting point with $x_1 < 0$ , a line search method that insists on staying interior cannot cross $x_1 = 0$ : doing so would require $x_2$ or $x_3$ to go negative. The iterates are driven toward the boundary of the interior region and stall at a point that is not a KKT point of the problem.

The remedy is the feasibility restoration phase, which we did not implement. Watch what our three strategies report --- and then watch what Ipopt, which does implement restoration, reports from the same starting point.

### Problem 6: the Wachter-Biegler counterexample

f6 = lambda x: x[0]


def c6(x):
    rhs = np.zeros(2)
    rhs[0] = x[0] ** 2 - x[1] - 1
    rhs[1] = x[0] - x[2] - 0.5
    return rhs


# Indices of variables with lower bound of zero
vb6 = [1, 2]

x0_6 = np.array([-2.0, 3.0, 1.0])

# Analytic solution
x_star_6 = np.array([1.0, 0.0, 0.5])
f_star_6 = 1.0

print("Analytic solution: x* =", x_star_6, " f(x*) =", f_star_6)
print()

results6 = compare_globalization(x0_6, f6, c6, vb6, reference=x_star_6)

Cross-check against Ipopt¶

Every answer above was checked against an analytic solution. Now check it a second way, with a production solver. Ipopt is an interior point method with a filter line search --- it is the mature version of what we just wrote --- so this is a meaningful comparison rather than a different algorithm agreeing by luck.

### Reference solutions from Ipopt

# This is the only cell in the assignment that uses Pyomo. Everything else is
# NumPy/SciPy, because the point of the assignment is to implement the
# algorithm, not to call one.

if "google.colab" in sys.modules:
    helper.install_idaes()
    helper.install_ipopt()

import pyomo.environ as pyo
from pyomo.opt import TerminationCondition


def solve_with_ipopt(build_model, name, expect_optimal=True):
    """Solve a Pyomo model with Ipopt and check the termination condition.

    Inputs:
        build_model - zero-argument function returning a ConcreteModel
        name - label for printing
        expect_optimal - if True, assert the solve succeeded

    Outputs:
        m - the solved model
        tc - the termination condition
    """

    m = build_model()
    results = pyo.SolverFactory("ipopt").solve(m, tee=False)
    tc = results.solver.termination_condition

    print("{0}: termination_condition = {1}".format(name, tc))

    if expect_optimal:
        assert tc == TerminationCondition.optimal, "Ipopt did not solve " + name

    return m, tc


def build4():
    m = pyo.ConcreteModel()
    m.x1 = pyo.Var(initialize=2.0)
    m.x2 = pyo.Var(initialize=2.0, bounds=(0, None))
    m.x3 = pyo.Var(initialize=1.0, bounds=(0, None))
    m.obj = pyo.Objective(expr=pyo.sqrt(1 + m.x1**2) + 0.5 * m.x1 + m.x2)
    m.c = pyo.Constraint(expr=m.x2 - m.x3 - 1 == 0)
    return m


def build5():
    m = pyo.ConcreteModel()
    m.x1 = pyo.Var(initialize=0.9)
    m.x2 = pyo.Var(initialize=0.9, bounds=(0, None))
    m.obj = pyo.Objective(expr=2 * m.x1 + m.x2)
    m.c = pyo.Constraint(expr=m.x1**2 + m.x2**2 - 1 == 0)
    return m


def build6(x1_init=-2.0):
    m = pyo.ConcreteModel()
    m.x1 = pyo.Var(initialize=x1_init)
    m.x2 = pyo.Var(initialize=3.0, bounds=(0, None))
    m.x3 = pyo.Var(initialize=1.0, bounds=(0, None))
    m.obj = pyo.Objective(expr=m.x1)
    m.c1 = pyo.Constraint(expr=m.x1**2 - m.x2 - 1 == 0)
    m.c2 = pyo.Constraint(expr=m.x1 - m.x3 - 0.5 == 0)
    return m


### Problems 4 and 5: Ipopt should agree with us and with the analytic solution

rows = []
for name, build, results, ours_f, f_star in [
    ("Problem 4", build4, results4, f4, f_star_4),
    ("Problem 5", build5, results5, f5, f_star_5),
]:
    m, tc = solve_with_ipopt(build, name)
    x_ours = np.asarray(results["filter"][0][-1], dtype=float)
    rows.append((name, ours_f(x_ours), pyo.value(m.obj), f_star))

print()
print("Problem   \tf(x), our filter\tf(x), Ipopt     \tf(x*), analytic")
print("-" * 84)
for name, ours, ipopt, exact in rows:
    print("{0} \t{1: 1.9e}\t{2: 1.9e}\t{3: 1.9e}".format(name, ours, ipopt, exact))

Problem 6: Ipopt fails here too¶

All three strategies failed on Problem 6, and it would be easy to blame our implementation. It is not our implementation. Ipopt is a mature, heavily tested filter line search interior point method with a feasibility restoration phase, and from the same starting point it reports infeasible and stops near the same non-KKT point our code stalls at.

That is why this problem is famous. Change the starting point and Ipopt finds the solution immediately --- the problem is not hard, the starting point is.

### Problem 6: Ipopt from the same starting point, and from a different one

m6a, tc6a = solve_with_ipopt(
    lambda: build6(x1_init=-2.0), "Problem 6 from x1 = -2", expect_optimal=False
)
print(
    "  Ipopt stops at      x  =",
    [round(pyo.value(z), 8) for z in (m6a.x1, m6a.x2, m6a.x3)],
)
print("  our filter stops at x  =", np.asarray(results6["filter"][0][-1], dtype=float))
print()

m6b, tc6b = solve_with_ipopt(lambda: build6(x1_init=2.0), "Problem 6 from x1 = +2")
print(
    "  Ipopt finds         x  =",
    [round(pyo.value(z), 8) for z in (m6b.x1, m6b.x2, m6b.x3)],
)
print("  analytic            x* =", x_star_6)
print()

assert tc6a == TerminationCondition.infeasible, (
    "Expected Ipopt to fail from x1 = -2; it did not. "
    "Check the Ipopt version before drawing conclusions from this cell."
)

print(
    "Both our code and Ipopt fail from x1 = -2 and succeed from x1 = +2.\n"
    "The failure is a property of the problem and the starting point,\n"
    "not of the implementation."
)

Discussion questions¶

Answer these using the numbers your notebook actually produced, not from memory.

Why does "none" diverge on Problem 4 but converge on Problems 1--3? Point to the specific structure of Problem 4 that breaks the full Newton step. Would a smaller starting value of $\mu_0$ have helped?
Problem 4 also runs a fourth variant that applies the step-to-the-boundary rule and then accepts $\alpha_{\max}^k$ with no acceptance test at all. It does not diverge, but it does not converge either. Explain what it is doing and why. What does this tell you about the division of labor between the step-to-the-boundary rule and the merit function or filter?
On Problem 5, the merit function and the filter disagree about the first step. Using the $(\theta, \varphi_{\mu_l})$ figure, explain the disagreement geometrically. Which acceptance region is larger, and is larger always better?
Count Newton iterations and trial points separately in the comparison table. A trial point costs one $f$ and one $c$ evaluation; a Newton iteration costs a Jacobian, a Hessian, and a linear solve. Which strategy is cheaper, and does the answer depend on the problem?
On Problem 6 every strategy fails, and so does Ipopt. Describe each failure mode --- what does "left the interior", "line search failed" and "restoration required" each mean here? Ipopt has a restoration phase and still fails from this starting point. What does that tell you about what restoration can and cannot fix?
The penalty update $\rho^k = \max\{\rho^{k-1}, \|v^k + d_v^k\|_\infty + \varrho\}$ is monotone --- $\rho$ never decreases. Print the $\rho$ history for Problem 5. What is the practical cost of an early overestimate of $\rho$ ?

Course policies and AI category¶

Tips and Tricks¶

Background¶

Problem Formulation¶

Reformulation Example¶

Primal Dual Optimality Conditions¶

Globalization¶

Step-to-the-boundary rule¶

Strategy 1: the ℓ1\ell_1ℓ1​ merit function¶

Strategy 2: the filter¶

Part 1: Basic Interior Point Method for Inequality and Equality Constraint NLPs¶

Pseudocode¶

Python Implementation¶

Adding globalization to the implementation¶

Test Problems¶

Problem 1: Convex¶

Problem 2: Convex¶

Problem 3: Nonconvex¶

Comparing Globalization Strategies¶

Problem 4: Where the full Newton step diverges¶

Problem 5: Where the merit function and the filter disagree¶

Problem 6: Where all three fail¶

Cross-check against Ipopt¶

Problem 6: Ipopt fails here too¶

Discussion questions¶

Strategy 1: the $\ell_1$ merit function¶