Convexity Revisited - Optimization for Decision Science

# This code cell installs packages on Colab

import sys

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

    helper.easy_install()
else:
    sys.path.insert(0, "../")
    import helper
helper.set_plotting_style()

import pyomo.environ as pyo

import random

# Seed the random number generator so this notebook is reproducible
# (Pyomo style guide, section 8).
random.seed(0)

Background¶

Reference: Beginning of Chapter 4 in Biegler (2010)

Canonical Nonlinear Program (NLP)¶

\begin{aligned} \min_{\mathbf{x}\in\mathbb R^n}\quad & f(\mathbf{x})\\ \text{s.t.}\quad & \mathbf{g}(\mathbf{x})\le\mathbf{0},\quad \mathbf{h}(\mathbf{x})=\mathbf{0}. \end{aligned}

(1)

The objective is scalar; $\mathbf{g}:\mathbb R^n\to\mathbb R^r$ and $\mathbf{h}:\mathbb R^n\to\mathbb R^m$ collect the constraints. For the derivative-based conditions below, assume continuous first and second derivatives near the point being tested. The feasible region is

\mathcal F=\{\mathbf{x}:\mathbf{g}(\mathbf{x})\le\mathbf{0},\ \mathbf{h}(\mathbf{x})=\mathbf{0}\}.

(2)

For a problem defined on a domain $X$ , intersect this set with $X$ . Each vector inequality is componentwise. Source: Biegler (2010), Eq. (4.1), p. 63.

Types of Constrained Optimal Solutions¶

All four definitions compare feasible points and require $\mathbf{x}^*\in\mathcal F$ .

Global: $f(\mathbf{x}^*)\le f(\mathbf{x})$ for every $\mathbf{x}\in\mathcal F$ .
Local: for some $\epsilon>0$ , the same inequality holds for every $\mathbf{x}\in\mathcal N_\epsilon(\mathbf{x}^*)\cap\mathcal F$ , where $\mathcal N_\epsilon=\{\mathbf{x}:\|\mathbf{x}-\mathbf{x}^*\|<\epsilon\}$ .
Strict local: there exists $\epsilon>0$ such that $f(\mathbf{x}^*)<f(\mathbf{x})$ for every $\mathbf{x}\in(\mathcal N_\epsilon(\mathbf{x}^*)\cap\mathcal F)\setminus\{\mathbf{x}^*\}$ .
Isolated local: $\mathbf{x}^*$ is local and there exists $\epsilon>0$ such that no distinct local minimizer belongs to $\mathcal N_\epsilon(\mathbf{x}^*)\cap\mathcal F$ .

Only the feasible part of a neighborhood supplies local comparisons.

The intersection is essential. For $f=x_1$ , $\mathcal F=\{x_1\ge0\}$ , the origin is a global minimum although its objective gradient is nonzero. It is not strict because feasible points on $x_1=0$ tie. For $f=x^2$ , $\mathcal F=[1,\infty)$ , the point 1 is strict, isolated and global. For $\mathcal F=\{0\}$ , the unique point is strict vacuously: there are no distinct feasible comparisons. For $f=0$ on $\mathbb R$ , every point is global/local but none is strict or isolated.

Three transfers of the minimizer definitions to different feasible sets.

The existing optimality reference gives the isolated-implies-strict argument and the strict-but-not-isolated counterexample; those proofs are not repeated here. Source: Biegler (2010), Definition 4.1, pp. 63–64, with the distinct-point qualification supplied.

Key Questions¶

As with unconstrained optimization, the following questions need to be considered:

If a solution $\mathbf{x}^*$ exists, is it a global solution in $\mathcal{F}$ or is it only a local solution?
What conditions characterize the optimal solutions?
Are there special problem classes of the NLP whose solutions have stronger properties and are easier to solve?
Are there efficient and reliable methods to solve the NLP?

Convexity for Constrained Optimization¶

Objective convexity and feasible-set convexity¶

A convex optimization problem has a convex feasible set and a convex objective on that set. Use a common convex domain $X$ when applying the function test:

f(\theta\mathbf{x}_a+(1-\theta)\mathbf{x}_b)\le \theta f(\mathbf{x}_a)+(1-\theta)f(\mathbf{x}_b),\qquad \mathbf{x}_a,\mathbf{x}_b\in X,\quad0<\theta<1.

(3)

Strict convexity makes this strict for distinct endpoints. We reserve $\theta$ for this mixture and $\alpha$ for a line-search length; the latter can exceed one.

The smooth bowl and upward jump are tested on explicitly marked domains.

For $(x^2-1)^2$ , the second derivative is $12x^2-4$ : it is positive on $[0.7,1.5]$ , but the larger $[-0.7,1.5]$ interval contains the comparison $f(0)=1>f(\pm0.7)=0.2601$ . The midpoint violates the chord bound. The second pictured function is $(\max\{0,1-x\})^2$ for $x\le2$ , then 1 for $x>2$ . On $[0,1.8]$ it is convex but not strict because it has a flat part. Across the jump, endpoints 2 and 2.5 have values 0 and 1, while their midpoint has value 1, greater than the chord value $1/2$ .

For the first interval of panel (b), the positive-part function $q(x)=\max\{0,1-x\}$ is convex and nonnegative. Thus, with $0<\theta<1$ ,

q(\theta a+(1-\theta)b)^2\le[\theta q(a)+(1-\theta)q(b)]^2 \le\theta q(a)^2+(1-\theta)q(b)^2.

(4)

The first inequality uses monotonicity of squaring on nonnegative values; the second is convexity of the scalar square. This proves convexity of the squared positive part without differentiating at the join. It is not strictly convex because it is constant on $[1,1.8]$ . For the quartic restricted to $[0.7,1.5]$ , $f''(x)\ge12(0.7)^2-4=1.88>0$ supplies the positive-curvature criterion from the Mathematics Primer.

A set is convex when it contains every segment between its points. Intersections preserve this property: endpoints in every set have their entire segment in every set. Unions need not preserve it; a midpoint between separated disks belongs to neither disk. The empty set is convex, but supplies no feasible solution.

A segment in an intersection, and a segment leaving a union.

Sources: Biegler (2010), Eq. (4.2) and Theorem 4.2, pp. 64–65; the explicit function examples are prepared versions of the handwritten domain concept test. The first-use definitions are in the Mathematics Primer.

Theorem 4.2: Convexity¶

Convex feasible-region certificate. On a common convex domain $X$ , suppose every scalar $g_i$ is convex and every equality is affine, $\mathbf h(\mathbf{x})=\mathbf A\mathbf{x}-\mathbf b$ . Then

\mathcal F=\{\mathbf{x}\in X:\mathbf g(\mathbf{x})\le0,\ \mathbf h(\mathbf{x})=0\}

(5)

is convex.

Take feasible endpoints and $\bar{\mathbf{x}}=\theta\mathbf{x}_a+(1-\theta)\mathbf{x}_b$ , with $0<\theta<1$ . Convexity of $X$ keeps $\bar{\mathbf{x}}$ in the domain. For each component,

g_i(\bar{\mathbf{x}})\le\theta g_i(\mathbf{x}_a)+(1-\theta)g_i(\mathbf{x}_b)\le0.

(6)

The affine equality gives

\mathbf h(\bar{\mathbf{x}})=\theta\mathbf h(\mathbf{x}_a)+(1-\theta)\mathbf h(\mathbf{x}_b)=0.

(7)

The constant is preserved because the weights sum to one. Thus the combination is feasible. In a contradiction proof, infeasibility means some $g_i(\bar{\mathbf{x}})>0$ or some equality fails; a vector expression $\mathbf g>0$ would incorrectly require every inequality component to be positive.

Constraint-by-constraint preservation of a feasible segment.

The equivalent set proof uses convex sublevel sets $\{g_i\le0\}$ and affine hyperplanes $\{h_j=0\}$ , intersecting all of them with $X$ .

Affine equalities are sufficient because they preserve both sides of an equality. Convexity alone gives only $h(\bar{\mathbf{x}})\le0$ . For $h=x_1^2+x_2^2-1$ , the endpoints $(-1,0),(1,0)$ satisfy $h=0$ , but their midpoint has $h=-1$ and is not on the circle.

Affine and curved zero-level sets give different midpoint tests.

Source: Biegler (2010), Theorem 4.2, pp. 64–65; the componentwise form follows the p.64 erratum. The affine extension follows directly from the displayed identity.

Theorem 4.3: Global Minimizers¶

Local-to-global certificate. If $\mathcal F$ is convex and $f$ is convex on it, every local minimizer is global. If $f$ is strictly convex, a minimizer, when it exists, is unique.

Suppose $\mathbf{x}_a$ is local but not global. Then some feasible point $\mathbf{x}_b$ has $f(\mathbf{x}_b)<f(\mathbf{x}_a)$ . It need not itself be a local minimizer. Let $\epsilon$ define the local neighborhood. For

\bar{\mathbf{x}}=\theta\mathbf{x}_a+(1-\theta)\mathbf{x}_b,

(8)

the segment is feasible and

\|\bar{\mathbf{x}}-\mathbf{x}_a\|=(1-\theta)\|\mathbf{x}_b-\mathbf{x}_a\|.

(9)

Choose $\theta$ sufficiently close to one to make this distance less than $\epsilon$ . Convexity gives

f(\bar{\mathbf{x}})\le f(\mathbf{x}_a)+(1-\theta) [f(\mathbf{x}_b)-f(\mathbf{x}_a)]<f(\mathbf{x}_a),

(10)

contradicting local minimality at a feasible neighbor.

The nearby comparison point remains in both the feasible set and local ball.

For uniqueness, suppose two distinct global minimizers tie at value $m$ . Their feasible midpoint has value strictly less than their average $m$ , a contradiction. Existence is separate: strictly convex $e^x$ on $\mathbb R$ has no minimizer.

These are applications of the existing first-use optimality lecture/reference argument to feasible comparisons. Sources: Biegler (2010), Theorem 4.3, p.65, correcting the proof’s two-local-minima assumption; Nocedal and Wright (2006), Theorem 2.5 proof, p.16.

More Illustrative Examples¶

For a nonempty bounded polytope, a linear objective attains an optimum at a vertex. An optimal face can contain other minimizers. This describes the bounded polygon pictures, not every LP: an LP can be infeasible, have an unbounded objective, or have a finite optimum on an unbounded polyhedron. A feasible polyhedron containing a line need not have vertices.

Strictly convex objectives can instead attain the unique minimum at a vertex, on an edge, or in the interior. Convexity certifies a value, not a universal boundary location.

The five convex contour cases used in the handout.

A nonconvex feasible region can supply separated local minima even when the objective is convex; a nonconvex objective can supply multiple local minima on a convex region. Prove a nonconvex-set claim by giving feasible endpoints and a segment point outside it.

In panel (1), the two local minima are at the inner-arc endpoints. In panel (2), they are at the centers of the three interior contour wells.

Sources: Biegler (2010), Figures 4.1–4.3, pp.66–67.

Circle Packing Example¶

Reference: Section 4.1 in Biegler (2010)

Motivating Question: Is this problem convex?

What is the smallest rectangle you can use to enclose three given circles? Reference: Example 4.4 in Biegler (2010).

Optimization Model and Pyomo Implementation¶

The following optimization model is given in Biegler (2010) and adapted to use set notation.

\begin{align*} \min_{\mathbf{x},\mathbf{y},A,B} \quad & 2(A + B) \\ \text{s.t.} \quad & A \geq 0, \quad B \geq 0 \\ & x_i, y_i \geq R_i, \quad x_i \leq B - R_i, \quad y_i \leq A - R_i, \quad \forall i \in \mathcal{C} \\ & (x_i - x_j)^2 + (y_i - y_j)^2 \geq (R_i + R_j)^2, \quad \forall i,j \in \{i \in \mathcal{C}, j \in \mathcal{C}: i < j\} \end{align*}

(11)

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches


def create_circle_model(circle_radii):
    """Create circle optimization model in Pyomo

    Arguments:
        circle_radii: dictionary with keys=circle name and value=radius (float)

    Returns:
        model: Pyomo model
    """

    # Create a concrete Pyomo model.
    model = pyo.ConcreteModel()

    # Initialize index for circles
    model.CIRCLES = pyo.Set(initialize=circle_radii.keys())

    # Create parameter
    model.R = pyo.Param(
        model.CIRCLES, domain=pyo.PositiveReals, initialize=circle_radii
    )

    # Create variables for box
    model.a = pyo.Var(domain=pyo.PositiveReals)
    model.b = pyo.Var(domain=pyo.PositiveReals)

    # Set objective
    model.obj = pyo.Objective(expr=2 * (model.a + model.b), sense=pyo.minimize)

    # Create variables for circle centers
    model.x = pyo.Var(model.CIRCLES, domain=pyo.PositiveReals)
    model.y = pyo.Var(model.CIRCLES, domain=pyo.PositiveReals)

    # "In the box" constraints
    @model.Constraint(model.CIRCLES)
    def left_x_con(m, c):
        return m.x[c] >= m.R[c]

    @model.Constraint(model.CIRCLES)
    def left_y_con(m, c):
        return m.y[c] >= m.R[c]

    @model.Constraint(model.CIRCLES)
    def right_x_con(m, c):
        return m.x[c] <= m.b - m.R[c]

    @model.Constraint(model.CIRCLES)
    def right_y_con(m, c):
        return m.y[c] <= m.a - m.R[c]

    # No overlap constraints
    @model.Constraint(model.CIRCLES, model.CIRCLES)
    def no_overlap_con(m, c1, c2):
        if c1 < c2:
            return (m.x[c1] - m.x[c2]) ** 2 + (m.y[c1] - m.y[c2]) ** 2 >= (
                m.R[c1] + m.R[c2]
            ) ** 2
        else:
            return pyo.Constraint.Skip

    return model


def initialize_circle_model(model, a_init=25, b_init=25):
    """Initialize the x and y coordinates using uniform distribution

    Arguments:
        model: Pyomo concrete model (modified in place)
        a_init: initial value for a (default=25)
        b_init: initial value for b (default=25)

    Returns:
        Nothing. But per Pyomo scoping rules, the input argument `model`
        can be modified in this function.

    """
    # Initialize
    model.a = a_init
    model.b = b_init

    for i in model.CIRCLES:
        # Adding circle radii ensures the remains in the >0, >0 quadrant
        model.x[i] = random.uniform(0, 10) + model.R[i]
        model.y[i] = random.uniform(0, 10) + model.R[i]

Next, we will create a dictionary containing the circle names and radii values.

# Create dictionary with circle data
circle_data = {"A": 10.0, "B": 5.0, "C": 3.0}
circle_data

{'A': 10.0, 'B': 5.0, 'C': 3.0}

# Access the keys
circle_data.keys()

dict_keys(['A', 'B', 'C'])

Now let’s create the model.

# Create model
model = create_circle_model(circle_data)

And let’s initialize the model.

# Initialize model
initialize_circle_model(model)
model.pprint()

1 Set Declarations
    CIRCLES : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    3 : {'A', 'B', 'C'}

1 Param Declarations
    R : Size=3, Index=CIRCLES, Domain=PositiveReals, Default=None, Mutable=False
        Key : Value
          A :  10.0
          B :   5.0
          C :   3.0

4 Var Declarations
    a : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :    25 :  None : False : False : PositiveReals
    b : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :    25 :  None : False : False : PositiveReals
    x : Size=3, Index=CIRCLES
        Key : Lower : Value              : Upper : Fixed : Stale : Domain
          A :     0 : 18.444218515250483 :  None : False : False : PositiveReals
          B :     0 :  9.205715808308451 :  None : False : False : PositiveReals
          C :     0 :  8.112747213686085 :  None : False : False : PositiveReals
    y : Size=3, Index=CIRCLES
        Key : Lower : Value              : Upper : Fixed : Stale : Domain
          A :     0 : 17.579544029403024 :  None : False : False : PositiveReals
          B :     0 :  7.589167502929634 :  None : False : False : PositiveReals
          C :     0 :  7.049341374504143 :  None : False : False : PositiveReals

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : 2*(a + b)

5 Constraint Declarations
    left_x_con : Size=3, Index=CIRCLES, Active=True
        Key : Lower : Body : Upper : Active
          A :  10.0 : x[A] :  +Inf :   True
          B :   5.0 : x[B] :  +Inf :   True
          C :   3.0 : x[C] :  +Inf :   True
    left_y_con : Size=3, Index=CIRCLES, Active=True
        Key : Lower : Body : Upper : Active
          A :  10.0 : y[A] :  +Inf :   True
          B :   5.0 : y[B] :  +Inf :   True
          C :   3.0 : y[C] :  +Inf :   True
    no_overlap_con : Size=3, Index=CIRCLES*CIRCLES, Active=True
        Key        : Lower : Body                                : Upper : Active
        ('A', 'B') : 225.0 : (x[A] - x[B])**2 + (y[A] - y[B])**2 :  +Inf :   True
        ('A', 'C') : 169.0 : (x[A] - x[C])**2 + (y[A] - y[C])**2 :  +Inf :   True
        ('B', 'C') :  64.0 : (x[B] - x[C])**2 + (y[B] - y[C])**2 :  +Inf :   True
    right_x_con : Size=3, Index=CIRCLES, Active=True
        Key : Lower : Body              : Upper : Active
          A :  -Inf : x[A] - (b - 10.0) :   0.0 :   True
          B :  -Inf :  x[B] - (b - 5.0) :   0.0 :   True
          C :  -Inf :  x[C] - (b - 3.0) :   0.0 :   True
    right_y_con : Size=3, Index=CIRCLES, Active=True
        Key : Lower : Body              : Upper : Active
          A :  -Inf : y[A] - (a - 10.0) :   0.0 :   True
          B :  -Inf :  y[B] - (a - 5.0) :   0.0 :   True
          C :  -Inf :  y[C] - (a - 3.0) :   0.0 :   True

12 Declarations: CIRCLES R a b obj x y left_x_con left_y_con right_x_con right_y_con no_overlap_con

Visualize Initial Point¶

Next, we’ll define a function to plot the solution (or initial point)

# Plot initial point


def plot_circles(m):
    """Plot circles using data in Pyomo model

    Arguments:
        m: Pyomo concrete model

    Returns:
        Nothing (but makes a figure)

    """

    # Create figure
    fig, ax = plt.subplots(1, figsize=(6, 6))

    # Adjust axes
    l = max(m.a.value, m.b.value) + 1
    ax.set_xlim(0, l)
    ax.set_ylim(0, l)

    # Draw box
    art = mpatches.Rectangle((0, 0), width=m.b.value, height=m.a.value, fill=False)
    ax.add_patch(art)

    # Draw circles and mark center
    for i in m.CIRCLES:
        art2 = mpatches.Circle(
            (m.x[i].value, m.y[i].value), radius=m.R[i], fill=True, alpha=0.25
        )
        ax.add_patch(art2)

        plt.scatter(m.x[i].value, m.y[i].value, color="black")

    # Show plot
    plt.show()


plot_circles(model)

Solve and Inspect the Solution¶

# Specify the solver
solver = pyo.SolverFactory("ipopt")

# Use the MA27 linear solver from HSL. The IDAES binaries
# (`idaes get-extensions`) are statically linked against HSL, so MA27 is
# available both here and on Colab. Pinning it keeps the iteration log in
# this notebook reproducible -- the default differs between Ipopt builds.
solver.options["linear_solver"] = "ma27"

# Solve the model
results = solver.solve(model, tee=True)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

Ipopt 3.13.2: linear_solver=ma27


******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:       30
Number of nonzeros in Lagrangian Hessian.............:       12

Total number of variables............................:        8
                     variables with only lower bounds:        8
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:       15
        inequality constraints with only lower bounds:        9
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        6

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e+02 6.25e+01 1.05e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.0046267e+02 4.59e+01 1.05e+01  -1.0 2.67e+01   0.0 1.52e-01 5.64e-02h  1
   2  1.0446269e+02 2.46e+01 2.54e+00  -1.0 5.43e+00  -0.5 5.62e-01 3.57e-01h  1
   3  1.1075999e+02 0.00e+00 4.40e+00  -1.0 4.83e+00  -1.0 6.87e-01 1.00e+00h  1
   4  1.0931790e+02 0.00e+00 4.22e+00  -1.0 6.13e+01  -1.4 3.27e-01 2.76e-02f  3
   5  1.0505265e+02 0.00e+00 8.83e-01  -1.0 2.84e+00  -1.0 1.00e+00 9.20e-01f  1
   6  1.0163361e+02 0.00e+00 1.15e+00  -1.0 1.96e+00  -1.5 4.90e-01 9.10e-01h  1
   7  9.9775211e+01 0.00e+00 1.90e-01  -1.0 2.12e+00  -1.1 1.00e+00 8.60e-01f  1
   8  9.8947266e+01 0.00e+00 8.86e-01  -1.0 1.23e+02    -  3.86e-01 5.70e-02f  1
   9  9.8937271e+01 0.00e+00 2.06e-02  -1.0 7.04e-01  -1.5 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  9.8325956e+01 0.00e+00 6.43e-03  -2.5 3.75e+00    -  9.93e-01 9.37e-01f  1
  11  9.8301250e+01 0.00e+00 2.66e-03  -2.5 4.81e+01    -  4.81e-01 1.00e+00h  1
  12  9.8301257e+01 0.00e+00 5.36e-04  -2.5 2.99e+01    -  1.00e+00 1.00e+00h  1
  13  9.8301256e+01 0.00e+00 2.96e-05  -2.5 2.59e+01    -  1.00e+00 1.00e+00h  1
  14  9.8301256e+01 0.00e+00 1.40e-06  -2.5 2.48e+00    -  1.00e+00 1.00e+00h  1
  15  9.8285159e+01 0.00e+00 2.69e-07  -3.8 3.80e-02    -  1.00e+00 1.00e+00h  1
  16  9.8284281e+01 0.00e+00 9.89e-10  -5.7 2.05e-03    -  1.00e+00 1.00e+00h  1
  17  9.8284270e+01 0.00e+00 9.89e-13  -8.6 1.62e-04    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 17

                                   (scaled)                 (unscaled)
Objective...............:   9.8284270438747257e+01    9.8284270438747257e+01
Dual infeasibility......:   9.8878042122038317e-13    9.8878042122038317e-13
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5165074041737747e-09    2.5165074041737747e-09
Overall NLP error.......:   2.5165074041737747e-09    2.5165074041737747e-09


Number of objective function evaluations             = 22
Number of objective gradient evaluations             = 18
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 22
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 18
Number of Lagrangian Hessian evaluations             = 17
Total CPU secs in IPOPT (w/o function evaluations)   =      0.002
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.

Next, we can inspect the solution. Because Pyomo is a Python extension, we can use Python (for loops, etc.) to programmatically inspect the solution.

# Print variable values
print("Name\tValue")
for c in model.component_data_objects(pyo.Var):
    print(c.name, "\t", pyo.value(c))

# Plot solution
plot_circles(model)

Name	Value
a 	 19.999999803189603
b 	 29.14213541618403
x[A] 	 19.14213551493119
x[B] 	 4.999999951252128
x[C] 	 5.517274811159322
y[A] 	 9.999999901937443
y[B] 	 4.999999953543125
y[C] 	 15.729315741741758

# Print constraints
for c in model.component_data_objects(pyo.Constraint):
    print(
        c.name,
        "\t",
        pyo.value(c.lower),
        "\t",
        pyo.value(c.body),
        "\t",
        pyo.value(c.upper),
    )

left_x_con[A] 	 10.0 	 19.14213551493119 	 None
left_x_con[B] 	 5.0 	 4.999999951252128 	 None
left_x_con[C] 	 3.0 	 5.517274811159322 	 None
left_y_con[A] 	 10.0 	 9.999999901937443 	 None
left_y_con[B] 	 5.0 	 4.999999953543125 	 None
left_y_con[C] 	 3.0 	 15.729315741741758 	 None
right_x_con[A] 	 None 	 9.87471615587765e-08 	 0.0
right_x_con[B] 	 None 	 -19.1421354649319 	 0.0
right_x_con[C] 	 None 	 -20.624860605024708 	 0.0
right_y_con[A] 	 None 	 9.874784012708915e-08 	 0.0
right_y_con[B] 	 None 	 -9.999999849646478 	 0.0
right_y_con[C] 	 None 	 -1.2706840614478452 	 0.0
no_overlap_con[A,B] 	 225.0 	 224.99999778541928 	 None
no_overlap_con[A,C] 	 169.0 	 218.46188918941948 	 None
no_overlap_con[B,C] 	 64.0 	 115.38579056358046 	 None

Reinitialize and Resolve¶

Reinitialize the model, plot the initial point, resolve, and plot the solution. Is there more than one solution?

# Initialize and print the model
initialize_circle_model(model)

# Plot initial point
plot_circles(model)

# Solve the model
results = solver.solve(model, tee=True)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

Ipopt 3.13.2: linear_solver=ma27


******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:       30
Number of nonzeros in Lagrangian Hessian.............:       12

Total number of variables............................:        8
                     variables with only lower bounds:        8
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:       15
        inequality constraints with only lower bounds:        9
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        6

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e+02 1.55e+02 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.0138375e+02 9.34e+01 7.73e-01  -1.0 1.19e+01    -  4.61e-01 3.31e-01h  1
   2  9.9526817e+01 4.07e+01 5.97e-01  -1.0 3.04e+01  -2.0 2.98e-01 1.77e-01h  1
   3  9.8395427e+01 2.47e+01 4.51e-01  -1.0 3.49e+00    -  9.69e-01 3.60e-01h  1
   4  9.8926255e+01 1.22e-01 9.29e-02  -1.0 3.37e+00    -  7.52e-01 8.39e-01h  1
   5  9.8741856e+01 5.57e-02 1.63e-01  -1.0 1.62e+02    -  7.54e-01 3.84e-01f  1
   6  9.8917368e+01 0.00e+00 3.14e-03  -1.0 3.30e+01    -  1.00e+00 1.00e+00h  1
   7  9.8900000e+01 0.00e+00 1.04e-03  -1.0 7.91e+00    -  1.00e+00 1.00e+00h  1
   8  9.8394513e+01 0.00e+00 4.20e-04  -1.7 1.76e+00    -  1.00e+00 1.00e+00h  1
   9  9.8284893e+01 0.00e+00 4.60e-05  -3.8 2.95e-01    -  1.00e+00 9.98e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  9.8284281e+01 0.00e+00 5.07e-08  -5.7 5.23e-02    -  1.00e+00 1.00e+00h  1
  11  9.8284270e+01 0.00e+00 3.77e-11  -8.6 4.78e-03    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 11

                                   (scaled)                 (unscaled)
Objective...............:   9.8284270438748450e+01    9.8284270438748450e+01
Dual infeasibility......:   3.7706208248288032e-11    3.7706208248288032e-11
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.7106270175648520e-09    2.7106270175648520e-09
Overall NLP error.......:   2.7106270175648520e-09    2.7106270175648520e-09


Number of objective function evaluations             = 12
Number of objective gradient evaluations             = 12
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 12
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 12
Number of Lagrangian Hessian evaluations             = 11
Total CPU secs in IPOPT (w/o function evaluations)   =      0.001
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.

# Plot solution
plot_circles(model)

Take Away Messages¶

Nonlinear programs may be nonconvex. For nonconvex problems, there often exist many local optima that are not also global optima.
Initialization is really important in optimization problems with nonlinear objectives or constraints!

Reference: representation tests and contour conclusions¶

The usual structural test checks a convex objective, convex scalar inequalities in $g_i\le0$ form, and affine equalities on a common convex domain. It is a sufficient certificate. A failed representation test does not alone prove the feasible set nonconvex.

For $\min x^2$ subject to $x^3\le0$ , the constraint function is not convex on $\mathbb R$ , since its second derivative $6x$ is negative for $x<0$ . Yet its feasible set is $(-\infty,0]$ , identical to the affine representation $x\le0$ . The problem has a convex objective and a convex feasible set. Similarly $x^3=0$ defines the singleton $\{0\}$ although that equality representation is nonlinear.

The two inequality representations give the same feasible half-line.

The circle-packing no-overlap set is outside a disk in relative-center coordinates. Opposite boundary points have an infeasible midpoint. The indexed model, two numerical arrangements and multistart experiment already have their home in Continuous Optimization: Nonlinear Programming; they are not copied into this reference. The existing numerical example above remains an introduction to the model.

Sources: Biegler (2010), Figures 4.1–4.3, pp.66–67, and Example 4.4, pp.67–68. The equivalent representations are directly verified by their signs/zero sets.