Polygon Clustering

Pack two rotatable polygons into the smallest axis-parallel rectangle without overlap, using a convex decomposition and nonlinear programming.

Notebook
Python / Knitro API Python / Pyomo Julia / JuMP

Introduction

Given two polygons (possibly non-convex) \(P_j\), \(j = 1, 2\), we look for an axis-parallel rectangle with the smallest area that contains the two polygons, without overlap between the polygons. Rotating the polygons is allowed.

This kind of problem arises in packing and nesting applications, for example when parts have to be cut out of a rectangular piece of stock material.

Here is an example of two input polygons and the clustering found by Knitro:

Problem description

Input

  • \(N\) polygons (here, \(N = 2\)), given by the coordinates of their vertices

Problem:

Choose, for each polygon:

  • a translation vector
  • a rotation angle

such that:

  • the polygons do not overlap
  • the polygons are contained in an axis-parallel rectangle

Objective:

Minimize the area of the containing rectangle.

Input data

A polygon is a ring of vertices with the first one repeated at the end, so the outline closes. centroid_of is the shoelace formula over that ring, and rotate turns a point about the origin given the sine and cosine of the angle rather than the angle itself, the form the model will carry a rotation. area is what the model minimizes, the width of the containing rectangle times its height.

A Packing holds that area together with the translation and the sine and cosine of each polygon, and that is all it takes to redraw the placement.

const Point = Tuple{Float64,Float64}

struct Vertex
    x::Float64
    y::Float64
end

struct ConvexPolygon
    vertices::Vector{Vertex}
    num_vertices::Int
end

function ConvexPolygon(ring::Vector{Point})
    ring[1] == ring[end] && (ring = ring[1:(end - 1)])
    return ConvexPolygon([Vertex(p[1], p[2]) for p in ring], length(ring))
end

mutable struct Polygon
    points::Vector{Point}
    centroid::Point
    parts::Vector{ConvexPolygon}
end

function Polygon(ring)
    points = Point[p for p in ring]
    return Polygon(points, centroid_of(points), ConvexPolygon[])
end

struct Packing
    area::Float64
    x::Vector{Float64}
    y::Vector{Float64}
    cosine::Vector{Float64}
    sine::Vector{Float64}
    width::Float64
    height::Float64
end

function centroid_of(points::Vector{Point})
    ring = points[1] == points[end] ? points[1:(end - 1)] : points
    twice, cx, cy = 0.0, 0.0, 0.0
    for i in eachindex(ring)
        x1, y1 = ring[i]
        x2, y2 = ring[mod1(i + 1, length(ring))]
        cross = x1 * y2 - x2 * y1
        twice += cross
        cx += (x1 + x2) * cross
        cy += (y1 + y2) * cross
    end
    return (cx / (3 * twice), cy / (3 * twice))
end

rotate(x, y, c, s) = (c * x + s * y, c * y - s * x)

area(placement) = placement.width * placement.height

Two polygons are packed below, both written out in full and both non-convex, the case the model has to handle.

polygons = [
    Polygon([
        (1, 0), (5, 0), (6, 3), (5, 7), (4, 4), (3, 8), (2, 3), (1, 6), (0, 3), (1, 0)
    ]),
    Polygon([
        (1, 0),
        (5, 0),
        (6, 4),
        (6, 7),
        (5, 4),
        (4, 7),
        (3, 2),
        (2, 6),
        (1, 3),
        (0, 7),
        (0, 3),
        (1, 0),
    ]),
]

A placement is drawn as the containing rectangle with the two outlines dropped into it, so what the picture shows is the gap the solver still has to close. The decomposition figure keeps the input outline on top of its convex parts, since the cuts are what the model sees and the outline is what the reader knows.

draw_polygons(polygons)
Polygon 1Polygon 2

Intersection of convex polygons

If two convex polygons are disjoint, then there is an edge \(e\) of one of them such that these polygons lie on the opposite sides of the line containing \(e\).

The position of a point relative to a line is what decides that. line_equation returns the coefficients \(a\), \(b\), \(c\) of the line through two points, and relative_position evaluates \(ax + by + c\). The sign says which side of the line the point falls on, so two points with opposite signs sit on opposite sides.

function line_equation(p1::Point, p2::Point)
    return (p2[2] - p1[2], p1[1] - p2[1], p2[1] * p1[2] - p1[1] * p2[2])
end

relative_position(a, b, c, point::Point) = a * point[1] + b * point[2] + c

Let’s illustrate this property:

triangle_1 = [(1.0, 2.0), (5.0, 5.0), (1.5, 6.0)]
triangle_2 = [(4.0, 3.4), (8.5, 0.6), (9.2, 3.2)]

draw_separation(triangle_1, triangle_2, (triangle_1[1], triangle_1[2]))
Two triangles and one separating line(AB)ABCDEF

Of the six edges these two triangles have, (AB) is the only one whose line separates them, with every vertex of one triangle staying in its own tint.

Running all six vertices through that line says the same thing in numbers. A and B are on it, the third vertex of the first triangle comes out negative, and all three of the second come out positive.

a, b, c = line_equation(triangle_1[1], triangle_1[2])

report_sides(a, b, c, [triangle_1, triangle_2])
vertex  triangle   ax + by + c  side       
-------------------------------------------
     A         1          0.00  on the line
     B         1          0.00  on the line
     C         1        -14.50  negative
     D         2          3.40  positive
     E         2         28.10  positive
     F         2         19.80  positive

Decomposition of the polygons into convex polygons

Since the intersection property only holds for convex polygons, in order to use it in our case, we need to decompose the input polygons into convex polygons first.

convex_parts cuts a ring into the fewest convex pieces it can. It takes a reflex vertex, tries every other vertex that one can see, and keeps the cut whose two halves need the fewest further cuts. The model works on those pieces rather than on the original ring.

Show the convex decomposition 
function triangle_area(a::Point, b::Point, c::Point)
    return (b[1] - a[1]) * (c[2] - a[2]) - (c[1] - a[1]) * (b[2] - a[2])
end

is_left_on(a, b, c) = triangle_area(a, b, c) >= 0
is_right(a, b, c) = triangle_area(a, b, c) < 0
is_right_on(a, b, c) = triangle_area(a, b, c) <= 0

square_distance(a::Point, b::Point) = (b[1] - a[1])^2 + (b[2] - a[2])^2

at(ring::Vector{Point}, i::Int) = ring[mod(i, length(ring)) + 1]

function line_crossing(a1::Point, a2::Point, b1::Point, b2::Point)
    ax, ay = a2[2] - a1[2], a1[1] - a2[1]
    ac = ax * a1[1] + ay * a1[2]
    bx, by = b2[2] - b1[2], b1[1] - b2[1]
    bc = bx * b1[1] + by * b1[2]
    det = ax * by - bx * ay
    det == 0 && return (0.0, 0.0)
    return ((by * ac - ay * bc) / det, (ax * bc - bx * ac) / det)
end

function is_reflex(ring::Vector{Point}, i::Int)
    return is_right(at(ring, i - 1), at(ring, i), at(ring, i + 1))
end

function can_see(ring::Vector{Point}, a::Int, b::Int)
    p, q = at(ring, a), at(ring, b)
    if is_left_on(at(ring, a + 1), p, q) && is_right_on(at(ring, a - 1), p, q)
        return false
    end
    reach = square_distance(p, q)
    for i in 0:(length(ring) - 1)
        (mod(i + 1, length(ring)) == a || i == a) && continue
        if is_left_on(p, q, at(ring, i + 1)) && is_right_on(p, q, at(ring, i))
            crossing = line_crossing(p, q, at(ring, i), at(ring, i + 1))
            square_distance(p, crossing) < reach && return false
        end
    end
    return true
end

function subring(ring::Vector{Point}, i::Int, j::Int)
    i < j && return ring[(i + 1):(j + 1)]
    return vcat(ring[1:(j + 1)], ring[(i + 1):end])
end

function cut_edges(ring::Vector{Point})
    best = Tuple{Point,Point}[]
    fewest = typemax(Int)
    for i in 0:(length(ring) - 1)
        is_reflex(ring, i) || continue
        for j in 0:(length(ring) - 1)
            can_see(ring, i, j) || continue
            cuts = vcat(cut_edges(subring(ring, i, j)), cut_edges(subring(ring, j, i)))
            if length(cuts) < fewest
                fewest = length(cuts)
                best = push!(cuts, (at(ring, i), at(ring, j)))
            end
        end
    end
    return best
end

function slice_at(ring::Vector{Point}, edge::Tuple{Point,Point})
    i = findfirst(==(edge[1]), ring)
    j = findfirst(==(edge[2]), ring)
    (i === nothing || j === nothing) && return nothing
    return (subring(ring, i - 1, j - 1), subring(ring, j - 1, i - 1))
end

function convex_parts(ring::Vector{Point})
    edges = cut_edges(ring)
    isempty(edges) && return [ring]
    parts = [ring]
    for edge in edges
        for (index, part) in enumerate(parts)
            pieces = slice_at(part, edge)
            if pieces !== nothing
                deleteat!(parts, index)
                append!(parts, [pieces[1], pieces[2]])
                break
            end
        end
    end
    return parts
end

function decompose!(polygon::Polygon)
    polygon.parts = ConvexPolygon.(convex_parts(polygon.points[1:(end - 1)]))
    return polygon
end
decompose!.(polygons)

draw_decompositions(polygons)
Polygon 1: 4 convex partsPolygon 2: 5 convex parts

Nonlinear model

Input

  • \(N\): number of polygons (here, \(N = 2\))
  • \(K^j\): number of convex parts of polygon \(P_j\)
  • \(L^j_k\): number of vertices of convex part \(k\) of polygon \(P_j\)
  • \(X^j_{k, l}\), \(Y^j_{k, l}\): x- and y-coordinates of vertex \(l\) of convex part \(k\) of polygon \(P_j\)

Variables

  • \(w, h \in \mathbb{R}^+\): width and height of the containing rectangle
  • \(x_j, y_j \in \mathbb{R}\), \(j = 1, \dots, N\): translation vector of polygon \(P_j\)
  • \(s_j, c_j \in [-1, 1]\), \(j = 1, \dots, N\): sine and cosine of the rotation angle of polygon \(P_j\)
  • \(x^j_{k, l}, y^j_{k, l} \in \mathbb{R}^+\): coordinates of vertex \(l\) of convex part \(k\) of polygon \(P_j\) after the rotation and the translation
  • \(d^{j_1, j_2}_{k_1, l_1, k_2} \in \{ 0, 1 \}\): 1 if the minimum relative position of the vertices of convex part \(k_2\) of polygon \(P_{j_2}\) with the line \((l_1, l_1 + 1)\) of convex part \(k_1\) of polygon \(P_{j_1}\) is non-negative

Objective: minimize the area of the containing rectangle

\[ \min w \cdot h \]

Constraints

  • Link between \(s_j\) and \(c_j\)

\[ \forall j = 1, \dots, N, \qquad c_j^2 + s_j^2 = 1 \]

  • Coordinates of the vertices

\[ \forall j = 1, \dots, N, \quad \forall k = 1, \dots, K^j, \quad \forall l = 1, \dots, L^j_k, \qquad x^j_{k, l} = x_j + X^j_{k, l} c_j + Y^j_{k, l} s_j \] \[ \forall j = 1, \dots, N, \quad \forall k = 1, \dots, K^j, \quad \forall l = 1, \dots, L^j_k, \qquad y^j_{k, l} = y_j + Y^j_{k, l} c_j - X^j_{k, l} s_j \]

  • Definition of \(d^{j_1, j_2}_{k_1, l_1, k_2}\)

\[ \forall j_1 = 1, \dots, N, \quad \forall j_2 = 1, \dots, N, j_2 \neq j_1, \quad \forall k_1 = 1, \dots, K^{j_1}, \quad \forall k_2 = 1, \dots, K^{j_2}, \quad \forall l_1 = 1, \dots, L^{j_1}_{k_1}, \quad \forall l_2 = 1, \dots, L^{j_2}_{k_2}, \]

\[ d^{j_1, j_2}_{k_1, l_1, k_2} \le \left( y^{j_1}_{k_1, l_1 + 1} - y^{j_1}_{k_1, l_1} \right) x^{j_2}_{k_2, l_2} + \left( x^{j_1}_{k_1, l_1} - x^{j_1}_{k_1, l_1 + 1} \right) y^{j_2}_{k_2, l_2} + x^{j_1}_{k_1, l_1 + 1} y^{j_1}_{k_1, l_1} - x^{j_1}_{k_1, l_1} y^{j_1}_{k_1, l_1 + 1} \]

  • For each pair of polygons, each pair of convex parts from each polygon must not intersect

\[ \forall j_1 = 1, \dots, N, \quad \forall j_2 = 1, \dots, N, j_2 \neq j_1, \quad \forall k_1 = 1, \dots, K^{j_1}, \quad \forall k_2 = 1, \dots, K^{j_2}, \]

\[ \max \left( \max_{l_1 = 1, \dots, L^{j_1}_{k_1}} d^{j_1, j_2}_{k_1, l_1, k_2}, \max_{l_2 = 1, \dots, L^{j_2}_{k_2}} d^{j_2, j_1}_{k_2, l_2, k_1} \right) \ge 0 \]

  • Link between \(w\) and \(h\), and \(x^j_{k, l}\) and \(y^j_{k, l}\)

\[ \forall j = 1, \dots, N, \quad \forall k = 1, \dots, K^j, \quad \forall l = 1, \dots, L^j_k, \qquad x^j_{k, l} \le w \] \[ \forall j = 1, \dots, N, \quad \forall k = 1, \dots, K^j, \quad \forall l = 1, \dots, L^j_k, \qquad y^j_{k, l} \le h \]

Solving the nonlinear model as it is wouldn’t give good results. To obtain good solutions quickly, we will:

  • first find a good initial point
  • then solve a simplified version of the model

Initial point

Starting from a good initial point is essential when solving packing problems with nonlinear programming.

generate_starts builds one start for every pair of rotations, eight of each polygon from 0° to 315° in 45° steps, so 64 in all. For each pair it moves the two polygons so their centroids coincide, slides the second one up until they barely overlap, and shifts the pair into the positive quadrant where the model expects it. separating_edge is what measures the overlap. It returns how deep the two convex parts still cut into each other, and the edge that keeps them apart best. That edge is also what the simplified model will fix.

const ROTATIONS = 0:45:315

struct Start
    angles::NTuple{2,Int}
    cosine::Vector{Float64}
    sine::Vector{Float64}
    x::Vector{Float64}
    y::Vector{Float64}
    width::Float64
    height::Float64
    separations::Vector{NTuple{5,Int}}
end

function Start(angles, x, y, width, height, edges)
    separations = [
        l1 != 0 ? (1, k1, l1, 2, k2) : (2, k2, l2, 1, k1) for (k1, k2, l1, l2) in edges
    ]
    cosine = [cos(deg2rad(angle)) for angle in angles]
    sine = [sin(deg2rad(angle)) for angle in angles]
    return Start(angles, cosine, sine, x, y, width, height, separations)
end

function separating_edge(part_1::Vector{Point}, part_2::Vector{Point})
    best, l1_best, l2_best = -Inf, 0, 0
    for l1 in eachindex(part_1)
        a, b, c = line_equation(part_1[l1], part_1[mod1(l1 + 1, length(part_1))])
        value = minimum(relative_position(a, b, c, v2) for v2 in part_2)
        if best < value
            best, l1_best, l2_best = value, l1, 0
        end
    end
    for l2 in eachindex(part_2)
        a, b, c = line_equation(part_2[l2], part_2[mod1(l2 + 1, length(part_2))])
        value = minimum(relative_position(a, b, c, v1) for v1 in part_1)
        if best < value
            best, l1_best, l2_best = value, 0, l2
        end
    end
    return -best, l1_best, l2_best
end

function worst_overlap(parts_1, parts_2)
    worst = 0.0
    edges = NTuple{4,Int}[]
    for (k1, part_1) in enumerate(parts_1), (k2, part_2) in enumerate(parts_2)
        depth, l1, l2 = separating_edge(part_1, part_2)
        worst = max(worst, depth)
        push!(edges, (k1, k2, l1, l2))
    end
    return worst, edges
end

function turned_parts(polygon::Polygon, theta)
    c, s = cos(deg2rad(theta)), sin(deg2rad(theta))
    xc, yc = rotate(polygon.centroid..., c, s)
    return [
        [
            (x - xc, y - yc) for
            (x, y) in (rotate(v.x, v.y, c, s) for v in part.vertices)
        ] for part in polygon.parts
    ]
end

function generate_starts(polygons::Vector{Polygon})
    starts = Start[]
    for theta1 in ROTATIONS
        parts_1 = turned_parts(polygons[1], theta1)
        for theta2 in ROTATIONS
            parts_2 = turned_parts(polygons[2], theta2)
            offset = 0.0
            local shifted, edges
            while true
                shifted = [[(x, y + offset) for (x, y) in part] for part in parts_2]
                worst, edges = worst_overlap(parts_1, shifted)
                worst < 1 && break
                offset += 1
            end

            points = vcat(Iterators.flatten(parts_1)..., Iterators.flatten(shifted)...)
            x_min = minimum(p[1] for p in points)
            y_min = minimum(p[2] for p in points)
            x_max = maximum(p[1] for p in points)
            y_max = maximum(p[2] for p in points)

            c1, s1 = cos(deg2rad(theta1)), sin(deg2rad(theta1))
            c2, s2 = cos(deg2rad(theta2)), sin(deg2rad(theta2))
            xc1, yc1 = rotate(polygons[1].centroid..., c1, s1)
            xc2, yc2 = rotate(polygons[2].centroid..., c2, s2)
            x = [-xc1 - x_min, -xc2 - x_min]
            y = [-yc1 - y_min, -yc2 + offset - y_min]
            width, height = x_max - x_min, y_max - y_min
            push!(starts, Start((theta1, theta2), x, y, width, height, edges))
        end
    end
    return sort!(starts; by=area)
end

The animation runs through all 64, first with the centroids on top of each other, then with the second polygon slid up clear of the first.

starts = generate_starts(polygons)
println("$(length(starts)) starts generated")

animate_starts(polygons, starts)
64 starts generated
1 / 64θ1 = 0° θ2 = 0°2 / 64θ1 = 0° θ2 = 45°3 / 64θ1 = 0° θ2 = 90°4 / 64θ1 = 0° θ2 = 135°5 / 64θ1 = 0° θ2 = 180°6 / 64θ1 = 0° θ2 = 225°7 / 64θ1 = 0° θ2 = 270°8 / 64θ1 = 0° θ2 = 315°9 / 64θ1 = 45° θ2 = 0°10 / 64θ1 = 45° θ2 = 45°11 / 64θ1 = 45° θ2 = 90°12 / 64θ1 = 45° θ2 = 135°13 / 64θ1 = 45° θ2 = 180°14 / 64θ1 = 45° θ2 = 225°15 / 64θ1 = 45° θ2 = 270°16 / 64θ1 = 45° θ2 = 315°17 / 64θ1 = 90° θ2 = 0°18 / 64θ1 = 90° θ2 = 45°19 / 64θ1 = 90° θ2 = 90°20 / 64θ1 = 90° θ2 = 135°21 / 64θ1 = 90° θ2 = 180°22 / 64θ1 = 90° θ2 = 225°23 / 64θ1 = 90° θ2 = 270°24 / 64θ1 = 90° θ2 = 315°25 / 64θ1 = 135° θ2 = 0°26 / 64θ1 = 135° θ2 = 45°27 / 64θ1 = 135° θ2 = 90°28 / 64θ1 = 135° θ2 = 135°29 / 64θ1 = 135° θ2 = 180°30 / 64θ1 = 135° θ2 = 225°31 / 64θ1 = 135° θ2 = 270°32 / 64θ1 = 135° θ2 = 315°33 / 64θ1 = 180° θ2 = 0°34 / 64θ1 = 180° θ2 = 45°35 / 64θ1 = 180° θ2 = 90°36 / 64θ1 = 180° θ2 = 135°37 / 64θ1 = 180° θ2 = 180°38 / 64θ1 = 180° θ2 = 225°39 / 64θ1 = 180° θ2 = 270°40 / 64θ1 = 180° θ2 = 315°41 / 64θ1 = 225° θ2 = 0°42 / 64θ1 = 225° θ2 = 45°43 / 64θ1 = 225° θ2 = 90°44 / 64θ1 = 225° θ2 = 135°45 / 64θ1 = 225° θ2 = 180°46 / 64θ1 = 225° θ2 = 225°47 / 64θ1 = 225° θ2 = 270°48 / 64θ1 = 225° θ2 = 315°49 / 64θ1 = 270° θ2 = 0°50 / 64θ1 = 270° θ2 = 45°51 / 64θ1 = 270° θ2 = 90°52 / 64θ1 = 270° θ2 = 135°53 / 64θ1 = 270° θ2 = 180°54 / 64θ1 = 270° θ2 = 225°55 / 64θ1 = 270° θ2 = 270°56 / 64θ1 = 270° θ2 = 315°57 / 64θ1 = 315° θ2 = 0°58 / 64θ1 = 315° θ2 = 45°59 / 64θ1 = 315° θ2 = 90°60 / 64θ1 = 315° θ2 = 135°61 / 64θ1 = 315° θ2 = 180°62 / 64θ1 = 315° θ2 = 225°63 / 64θ1 = 315° θ2 = 270°64 / 64θ1 = 315° θ2 = 315°
Centroids coincide
area = 6.04 × 15.15 = 91.55area = 8.49 × 16.08 = 136.42area = 7.35 × 13.98 = 102.71area = 8.49 × 15.32 = 130.01area = 6.17 × 10.67 = 65.85smallest of the 64area = 8.49 × 16.23 = 137.71area = 7.13 × 13.85 = 98.78area = 8.49 × 15.99 = 135.64area = 7.78 × 14.64 = 113.86area = 8.62 × 15.56 = 134.20area = 7.79 × 13.46 = 104.87area = 8.72 × 15.81 = 137.76area = 7.78 × 14.15 = 110.10area = 9.30 × 12.71 = 118.31area = 8.47 × 13.33 = 112.94area = 9.40 × 15.47 = 145.38area = 8.15 × 13.13 = 107.11area = 8.50 × 14.06 = 119.51area = 8.00 × 11.96 = 95.67area = 8.49 × 14.30 = 121.37area = 8.02 × 12.65 = 101.51area = 10.16 × 14.21 = 144.44area = 9.33 × 11.83 = 110.36area = 10.26 × 14.97 = 153.51area = 7.07 × 15.47 = 109.40area = 8.49 × 16.40 = 139.13area = 7.64 × 13.30 = 101.56area = 8.56 × 15.64 = 133.94area = 7.07 × 14.99 = 105.98area = 8.75 × 16.55 = 144.78area = 7.92 × 14.16 = 112.13area = 8.84 × 16.30 = 144.15area = 6.17 × 15.33 = 94.63area = 8.49 × 16.26 = 137.94area = 7.13 × 14.15 = 100.99area = 8.49 × 16.50 = 140.01area = 6.04 × 14.85 = 89.71area = 8.49 × 16.41 = 139.22area = 7.35 × 14.02 = 103.06area = 8.49 × 17.16 = 145.64area = 7.78 × 13.92 = 108.24area = 9.30 × 14.84 = 138.09area = 8.47 × 12.74 = 107.93area = 9.40 × 15.09 = 141.76area = 7.78 × 13.43 = 104.48area = 8.62 × 14.99 = 129.28area = 7.79 × 12.61 = 98.22area = 8.72 × 15.75 = 137.25area = 8.02 × 13.35 = 107.11area = 10.16 × 14.27 = 145.08area = 9.33 × 12.17 = 113.58area = 10.26 × 14.52 = 148.90area = 8.15 × 12.87 = 104.91area = 8.50 × 14.43 = 122.62area = 8.00 × 12.04 = 96.33area = 8.49 × 14.18 = 120.34area = 7.07 × 14.79 = 104.58area = 8.75 × 15.72 = 137.49area = 7.92 × 13.61 = 107.77area = 8.84 × 14.96 = 132.26area = 7.07 × 14.31 = 101.16area = 8.49 × 15.87 = 134.63area = 7.64 × 14.48 = 110.63area = 8.56 × 15.62 = 133.79
Polygon 2 slides up

Laid out as a grid, the same 64 starts show where the room goes. The row is the turn applied to polygon 1 and the column the turn applied to polygon 2, every cell is drawn to the same scale so the rectangles compare by eye, and the number under each one is its area.

draw_start_grid(polygons, starts)
Every start, by rotation pairPolygon 2 turned byPolygon 1 turned by0°45°90°135°180°225°270°315°0°45°90°135°180°225°270°315°θ1 = 0°, θ2 = 0°, 6.04 × 15.15, area 91.5591.5θ1 = 0°, θ2 = 45°, 8.49 × 16.08, area 136.42136.4θ1 = 0°, θ2 = 90°, 7.35 × 13.98, area 102.71102.7θ1 = 0°, θ2 = 135°, 8.49 × 15.32, area 130.01130.0θ1 = 0°, θ2 = 180°, 6.17 × 10.67, area 65.8565.9 ✓θ1 = 0°, θ2 = 225°, 8.49 × 16.23, area 137.71137.7θ1 = 0°, θ2 = 270°, 7.13 × 13.85, area 98.7898.8θ1 = 0°, θ2 = 315°, 8.49 × 15.99, area 135.64135.6θ1 = 45°, θ2 = 0°, 7.78 × 14.64, area 113.86113.9θ1 = 45°, θ2 = 45°, 8.62 × 15.56, area 134.20134.2θ1 = 45°, θ2 = 90°, 7.79 × 13.46, area 104.87104.9θ1 = 45°, θ2 = 135°, 8.72 × 15.81, area 137.76137.8θ1 = 45°, θ2 = 180°, 7.78 × 14.15, area 110.10110.1θ1 = 45°, θ2 = 225°, 9.30 × 12.71, area 118.31118.3θ1 = 45°, θ2 = 270°, 8.47 × 13.33, area 112.94112.9θ1 = 45°, θ2 = 315°, 9.40 × 15.47, area 145.38145.4θ1 = 90°, θ2 = 0°, 8.15 × 13.13, area 107.11107.1θ1 = 90°, θ2 = 45°, 8.50 × 14.06, area 119.51119.5θ1 = 90°, θ2 = 90°, 8.00 × 11.96, area 95.6795.7θ1 = 90°, θ2 = 135°, 8.49 × 14.30, area 121.37121.4θ1 = 90°, θ2 = 180°, 8.02 × 12.65, area 101.51101.5θ1 = 90°, θ2 = 225°, 10.16 × 14.21, area 144.44144.4θ1 = 90°, θ2 = 270°, 9.33 × 11.83, area 110.36110.4θ1 = 90°, θ2 = 315°, 10.26 × 14.97, area 153.51153.5θ1 = 135°, θ2 = 0°, 7.07 × 15.47, area 109.40109.4θ1 = 135°, θ2 = 45°, 8.49 × 16.40, area 139.13139.1θ1 = 135°, θ2 = 90°, 7.64 × 13.30, area 101.56101.6θ1 = 135°, θ2 = 135°, 8.56 × 15.64, area 133.94133.9θ1 = 135°, θ2 = 180°, 7.07 × 14.99, area 105.98106.0θ1 = 135°, θ2 = 225°, 8.75 × 16.55, area 144.78144.8θ1 = 135°, θ2 = 270°, 7.92 × 14.16, area 112.13112.1θ1 = 135°, θ2 = 315°, 8.84 × 16.30, area 144.15144.2θ1 = 180°, θ2 = 0°, 6.17 × 15.33, area 94.6394.6θ1 = 180°, θ2 = 45°, 8.49 × 16.26, area 137.94137.9θ1 = 180°, θ2 = 90°, 7.13 × 14.15, area 100.99101.0θ1 = 180°, θ2 = 135°, 8.49 × 16.50, area 140.01140.0θ1 = 180°, θ2 = 180°, 6.04 × 14.85, area 89.7189.7θ1 = 180°, θ2 = 225°, 8.49 × 16.41, area 139.22139.2θ1 = 180°, θ2 = 270°, 7.35 × 14.02, area 103.06103.1θ1 = 180°, θ2 = 315°, 8.49 × 17.16, area 145.64145.6θ1 = 225°, θ2 = 0°, 7.78 × 13.92, area 108.24108.2θ1 = 225°, θ2 = 45°, 9.30 × 14.84, area 138.09138.1θ1 = 225°, θ2 = 90°, 8.47 × 12.74, area 107.93107.9θ1 = 225°, θ2 = 135°, 9.40 × 15.09, area 141.76141.8θ1 = 225°, θ2 = 180°, 7.78 × 13.43, area 104.48104.5θ1 = 225°, θ2 = 225°, 8.62 × 14.99, area 129.28129.3θ1 = 225°, θ2 = 270°, 7.79 × 12.61, area 98.2298.2θ1 = 225°, θ2 = 315°, 8.72 × 15.75, area 137.25137.3θ1 = 270°, θ2 = 0°, 8.02 × 13.35, area 107.11107.1θ1 = 270°, θ2 = 45°, 10.16 × 14.27, area 145.08145.1θ1 = 270°, θ2 = 90°, 9.33 × 12.17, area 113.58113.6θ1 = 270°, θ2 = 135°, 10.26 × 14.52, area 148.90148.9θ1 = 270°, θ2 = 180°, 8.15 × 12.87, area 104.91104.9θ1 = 270°, θ2 = 225°, 8.50 × 14.43, area 122.62122.6θ1 = 270°, θ2 = 270°, 8.00 × 12.04, area 96.3396.3θ1 = 270°, θ2 = 315°, 8.49 × 14.18, area 120.34120.3θ1 = 315°, θ2 = 0°, 7.07 × 14.79, area 104.58104.6θ1 = 315°, θ2 = 45°, 8.75 × 15.72, area 137.49137.5θ1 = 315°, θ2 = 90°, 7.92 × 13.61, area 107.77107.8θ1 = 315°, θ2 = 135°, 8.84 × 14.96, area 132.26132.3θ1 = 315°, θ2 = 180°, 7.07 × 14.31, area 101.16101.2θ1 = 315°, θ2 = 225°, 8.49 × 15.87, area 134.63134.6θ1 = 315°, θ2 = 270°, 7.64 × 14.48, area 110.63110.6θ1 = 315°, θ2 = 315°, 8.56 × 15.62, area 133.79133.8

Rotations that interlock the two outlines need the least room. Here is the smallest of the 64, the start the solver gets:

start = starts[1]
draw_placement(polygons, start, "Smallest start")
Smallest start: 6.17 × 10.67, area 65.85

That start already looks close, but it is not a solution. The polygons still overlap slightly, and they do not sit snugly against each other. Closing that gap is what the nonlinear program is for.

Simplified nonlinear model

To simplify the model, for each pair of convex polygons from different polygons, we will fix in advance which edge from these convex polygons is responsible for the non-intersection of the convex polygons. The separating edges are chosen when generating the initial point: they are the edges that separate the two convex polygons the most in the initial placement.

That is, we replace these constraints

\[ \forall j_1 = 1, \dots, N, \quad \forall j_2 = 1, \dots, N, j_2 \neq j_1, \quad \forall k_1 = 1, \dots, K^{j_1}, \quad \forall k_2 = 1, \dots, K^{j_2}, \quad \forall l_1 = 1, \dots, L^{j_1}_{k_1}, \quad \forall l_2 = 1, \dots, L^{j_2}_{k_2}, \]

\[ d^{j_1, j_2}_{k_1, l_1, k_2} \le \left( y^{j_1}_{k_1, l_1 + 1} - y^{j_1}_{k_1, l_1} \right) x^{j_2}_{k_2, l_2} + \left( x^{j_1}_{k_1, l_1} - x^{j_1}_{k_1, l_1 + 1} \right) y^{j_2}_{k_2, l_2} + x^{j_1}_{k_1, l_1 + 1} y^{j_1}_{k_1, l_1} - x^{j_1}_{k_1, l_1} y^{j_1}_{k_1, l_1 + 1} \]

and

\[ \forall j_1 = 1, \dots, N, \quad \forall j_2 = 1, \dots, N, j_2 \neq j_1, \quad \forall k_1 = 1, \dots, K^{j_1}, \quad \forall k_2 = 1, \dots, K^{j_2}, \]

\[ \max \left( \max_{l_1 = 1, \dots, L^{j_1}_{k_1}} d^{j_1, j_2}_{k_1, l_1, k_2}, \max_{l_2 = 1, \dots, L^{j_2}_{k_2}} d^{j_2, j_1}_{k_2, l_2, k_1} \right) \ge 0 \]

by the following constraints. If the edge chosen to separate convex part \(k_1\) of polygon \(P_{j_1}\) and convex part \(k_2\) of polygon \(P_{j_2}\) is edge \((l_1, l_1 + 1)\) from convex part \(k_1\) of polygon \(P_{j_1}\):

\[ \forall l_2 = 1, \dots, L^{j_2}_{k_2}, \qquad \left( y^{j_1}_{k_1, l_1 + 1} - y^{j_1}_{k_1, l_1} \right) x^{j_2}_{k_2, l_2} + \left( x^{j_1}_{k_1, l_1} - x^{j_1}_{k_1, l_1 + 1} \right) y^{j_2}_{k_2, l_2} + x^{j_1}_{k_1, l_1 + 1} y^{j_1}_{k_1, l_1} - x^{j_1}_{k_1, l_1} y^{j_1}_{k_1, l_1 + 1} \ge 0 \]

Otherwise, if the edge chosen to separate convex part \(k_1\) of polygon \(P_{j_1}\) and convex part \(k_2\) of polygon \(P_{j_2}\) is edge \((l_2, l_2 + 1)\) from convex part \(k_2\) of polygon \(P_{j_2}\):

\[ \forall l_1 = 1, \dots, L^{j_1}_{k_1}, \qquad \left( y^{j_2}_{k_2, l_2 + 1} - y^{j_2}_{k_2, l_2} \right) x^{j_1}_{k_1, l_1} + \left( x^{j_2}_{k_2, l_2} - x^{j_2}_{k_2, l_2 + 1} \right) y^{j_1}_{k_1, l_1} + x^{j_2}_{k_2, l_2 + 1} y^{j_2}_{k_2, l_2} - x^{j_2}_{k_2, l_2} y^{j_2}_{k_2, l_2 + 1} \ge 0 \]

Variables \(d^{j_1, j_2}_{k_1, l_1, k_2}\) are not necessary anymore: the simplified model is a pure nonlinear program (NLP).

Model implementation

using JuMP
using KNITRO

minimize_area(polygons, start) builds the simplified model, solves it with Knitro, and returns a Packing. Every expression in it is quadratic, from the area \(w \cdot h\) to the circle \(c_j^2 + s_j^2 = 1\) that keeps the rotation a rotation, and the separation constraints.

The start is only nearly feasible, so the starting vertex coordinates are clamped to the bounds the model gives them.

function minimize_area(polygons::Vector{Polygon}, start::Start)
    num_polygons = length(polygons)

    vertices = [
        (j, k, i) for (j, polygon) in enumerate(polygons) for
        (k, part) in enumerate(polygon.parts) for i in 1:part.num_vertices
    ]

    model = Model(KNITRO.Optimizer)

    @variable(model, w >= 0, start = start.width)
    @variable(model, h >= 0, start = start.height)
    @variable(model, x[j in 1:num_polygons], start = start.x[j])
    @variable(model, y[j in 1:num_polygons], start = start.y[j])
    @variable(model, -1 <= c[j in 1:num_polygons] <= 1, start = start.cosine[j])
    @variable(model, -1 <= s[j in 1:num_polygons] <= 1, start = start.sine[j])
    @variable(model, xv[v in vertices] >= 0)
    @variable(model, yv[v in vertices] >= 0)
    for (j, k, i) in vertices
        v = polygons[j].parts[k].vertices[i]
        v_x, v_y = rotate(v.x, v.y, start.cosine[j], start.sine[j])
        set_start_value(xv[(j, k, i)], max(0.0, start.x[j] + v_x))
        set_start_value(yv[(j, k, i)], max(0.0, start.y[j] + v_y))
    end

    @objective(model, Min, w * h)

    for j in 1:num_polygons
        @constraint(model, c[j]^2 + s[j]^2 == 1)
    end

    for (j, k, i) in vertices
        v = polygons[j].parts[k].vertices[i]
        @constraint(model, xv[(j, k, i)] == x[j] + v.x * c[j] + v.y * s[j])
        @constraint(model, yv[(j, k, i)] == y[j] + v.y * c[j] - v.x * s[j])
        @constraint(model, xv[(j, k, i)] <= w)
        @constraint(model, yv[(j, k, i)] <= h)
    end

    for (j1, k1, l1, j2, k2) in start.separations
        p1 = (j1, k1, l1)
        p2 = (j1, k1, mod1(l1 + 1, polygons[j1].parts[k1].num_vertices))
        for l2 in 1:(polygons[j2].parts[k2].num_vertices)
            v = (j2, k2, l2)
            expr = (yv[p2] - yv[p1]) * xv[v]
            expr += (xv[p1] - xv[p2]) * yv[v]
            expr += xv[p2] * yv[p1] - xv[p1] * yv[p2]
            @constraint(model, expr >= 0)
        end
    end

    optimize!(model)

    smallest_area = objective_value(model)
    solved = [value.(vector) for vector in (x, y, c, s)]
    return Packing(smallest_area, solved..., value(w), value(h))
end

Resolution

Let’s solve the simplified model from the best initial point.

packing = minimize_area(polygons, start);
=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  216 |                           216
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        78         0 |         2        78         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       472       345 |         2       472       345
  total:              2       852       345 |         2       852       345
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    6.585371e+01   1.94e-01
       4    6.000000e+01   3.12e-12   9.21e-09   4.07e-05      0.02

EXIT: Locally optimal solution found.

HINT: Knitro spent   2.7% of solution time (0.000486 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   6.00000000101676e+01
Final feasibility error (abs / rel) =   3.12e-12 / 3.12e-12
Final optimality error  (abs / rel) =   9.21e-09 / 9.21e-10
# of iterations                     =          4 
# of CG iterations                  =          0 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.01808 (     0.013 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================
Show the full outputHide the full output

Output visualization

Let’s visualize the solution found by the solver.

report_packing(packing)
Rectangle 6.000 × 10.000, area 60.000

polygon         x         y      cos      sin     turn
------------------------------------------------------
      1     0.000     0.000    1.000   -0.000     0.0°
      2     6.000    10.000   -1.000   -0.000   180.0°
draw_placement(polygons, packing, "Knitro solution")
Knitro solution: 6.00 × 10.00, area 60.00

Resolution from each initial point

The simplified model is a nonconvex NLP, so the solution found by the solver depends on the initial point. Each solve takes a fraction of a second, so all 64 starts can be run and compared.

packings = [minimize_area(polygons, candidate) for candidate in starts];
Show the full outputHide the full output
=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  216 |                           216
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        78         0 |         2        78         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       472       345 |         2       472       345
  total:              2       852       345 |         2       852       345
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    6.585371e+01   1.94e-01
       4    6.000000e+01   3.12e-12   9.21e-09   4.07e-05      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   7.0% of solution time (0.000369 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   6.00000000101676e+01
Final feasibility error (abs / rel) =   3.12e-12 / 3.12e-12
Final optimality error  (abs / rel) =   9.21e-09 / 9.21e-10
# of iterations                     =          4 
# of CG iterations                  =          0 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00533 (     0.005 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  214 |                           214
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        76         0 |         2        76         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       460       333 |         2       460       333
  total:              2       840       333 |         2       840       333
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    8.970511e+01   6.10e-01
       6    8.864894e+01   6.43e-10   7.36e-08   6.35e-05      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   2.7% of solution time (0.000221 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   8.86489397746278e+01
Final feasibility error (abs / rel) =   6.43e-10 / 6.43e-10
Final optimality error  (abs / rel) =   7.36e-08 / 4.98e-09
# of iterations                     =          6 
# of CG iterations                  =          1 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00809 (     0.008 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  217 |                           217
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        79         0 |         2        79         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       478       339 |         2       478       339
  total:              2       858       339 |         2       858       339
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    9.154839e+01   1.78e-15
       3    9.000000e+01   8.78e-10   5.71e-06   9.91e-05      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   4.4% of solution time (0.000223 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.00000004808873e+01
Final feasibility error (abs / rel) =   8.78e-10 / 8.78e-10
Final optimality error  (abs / rel) =   5.71e-06 / 3.81e-07
# of iterations                     =          3 
# of CG iterations                  =          0 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00511 (     0.005 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  217 |                           217
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        79         0 |         2        79         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       478       337 |         2       478       337
  total:              2       858       337 |         2       858       337
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    9.463135e+01   1.78e-15
       4    9.000000e+01   2.00e-11   4.01e-07   3.72e-05      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   3.8% of solution time (0.000234 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.00000000750495e+01
Final feasibility error (abs / rel) =   2.00e-11 / 2.00e-11
Final optimality error  (abs / rel) =   4.01e-07 / 2.68e-08
# of iterations                     =          4 
# of CG iterations                  =          0 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00613 (     0.006 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  211 |                           211
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        73         0 |         2        73         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       442       315 |         2       442       315
  total:              2       822       315 |         2       822       315
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    9.566573e+01   4.44e-16
      10    9.200038e+01   3.62e-07   7.16e-01   4.08e-03      0.01
      12    9.200014e+01   4.85e-08   1.37e-06   3.43e-04      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   1.6% of solution time (0.000233 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.20001376874382e+01
Final feasibility error (abs / rel) =   4.85e-08 / 4.85e-08
Final optimality error  (abs / rel) =   1.37e-06 / 1.19e-07
# of iterations                     =         12 
# of CG iterations                  =          3 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.01436 (     0.014 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  213 |                           213
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        75         0 |         2        75         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       454       323 |         2       454       323
  total:              2       834       323 |         2       834       323
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    9.633427e+01   1.78e-15
      10    9.333810e+01   9.28e-07   3.21e+00   1.14e-02      0.01
      20    9.330939e+01   8.00e-10   4.32e+00   8.42e-06      0.16
      30    9.650491e+01   1.06e-01   2.46e+00   1.54e+00      0.17
      40    9.330955e+01   1.74e-08   6.59e-04   8.57e-04      0.17
      42    9.330938e+01   7.78e-11   1.77e-10   1.89e-07      0.18

EXIT: Locally optimal solution found.

Final Statistics
----------------
Final objective value               =   9.33093779780084e+01
Final feasibility error (abs / rel) =   7.78e-11 / 7.78e-11
Final optimality error  (abs / rel) =   1.77e-10 / 1.52e-11
# of iterations                     =         42 
# of CG iterations                  =          7 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.17514 (     0.032 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  212 |                           212
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        74         0 |         2        74         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       448       315 |         2       448       315
  total:              2       828       315 |         2       828       315
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    9.822454e+01   4.24e-01
      10    9.832886e+01   4.21e-07   1.74e-04   3.71e-03      0.01
      20    9.813334e+01   2.18e-07   4.56e-02   8.43e-03      0.02
      21    9.813333e+01   1.69e-09   1.01e-05   6.19e-06      0.02

EXIT: Locally optimal solution found.

HINT: Knitro spent   1.0% of solution time (0.000234 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.81333333529302e+01
Final feasibility error (abs / rel) =   1.69e-09 / 1.69e-09
Final optimality error  (abs / rel) =   1.01e-05 / 8.20e-07
# of iterations                     =         21 
# of CG iterations                  =         15 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.02303 (     0.023 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  213 |                           213
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        75         0 |         2        75         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       454       321 |         2       454       321
  total:              2       834       321 |         2       834       321
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    9.878228e+01   3.29e-01
       9    9.318974e+01   3.28e-10   4.51e-07   2.46e-05      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   1.9% of solution time (0.000230 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.31897436123697e+01
Final feasibility error (abs / rel) =   3.28e-10 / 3.28e-10
Final optimality error  (abs / rel) =   4.51e-07 / 3.39e-08
# of iterations                     =          9 
# of CG iterations                  =          3 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.01193 (     0.012 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  217 |                           217
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        79         0 |         2        79         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       478       337 |         2       478       337
  total:              2       858       337 |         2       858       337
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.009890e+02   4.44e-16
       5    9.569454e+01   2.67e-10   9.88e-08   3.79e-04      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   3.7% of solution time (0.000226 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.56945403904340e+01
Final feasibility error (abs / rel) =   2.67e-10 / 2.67e-10
Final optimality error  (abs / rel) =   9.88e-08 / 7.23e-09
# of iterations                     =          5 
# of CG iterations                  =          0 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00607 (     0.006 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  214 |                           214
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        76         0 |         2        76         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       460       321 |         2       460       321
  total:              2       840       321 |         2       840       321
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.011610e+02   3.56e-01
       5    1.004680e+02   1.63e-11   3.94e-09   1.71e-05      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   3.2% of solution time (0.000223 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   1.00468037437097e+02
Final feasibility error (abs / rel) =   1.63e-11 / 1.63e-11
Final optimality error  (abs / rel) =   3.94e-09 / 2.77e-10
# of iterations                     =          5 
# of CG iterations                  =          3 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00700 (     0.007 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  213 |                           213
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        75         0 |         2        75         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       454       321 |         2       454       321
  total:              2       834       321 |         2       834       321
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.015109e+02   2.30e-02
       3    1.013333e+02   1.18e-10   2.71e-06   1.03e-03      0.00

EXIT: Locally optimal solution found.

HINT: Knitro spent   4.5% of solution time (0.000223 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   1.01333337801904e+02
Final feasibility error (abs / rel) =   1.18e-10 / 1.18e-10
Final optimality error  (abs / rel) =   2.71e-06 / 2.14e-07
# of iterations                     =          3 
# of CG iterations                  =          0 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00495 (     0.005 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  211 |                           211
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        73         0 |         2        73         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       442       313 |         2       442       313
  total:              2       822       313 |         2       822       313
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.015563e+02   1.21e-01
       6    9.522971e+01   1.07e-11   1.60e-06   4.80e-06      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   3.2% of solution time (0.000218 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.52297077513258e+01
Final feasibility error (abs / rel) =   1.07e-11 / 1.07e-11
Final optimality error  (abs / rel) =   1.60e-06 / 1.19e-07
# of iterations                     =          6 
# of CG iterations                  =          0 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00677 (     0.007 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  213 |                           213
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        75         0 |         2        75         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       454       315 |         2       454       315
  total:              2       834       315 |         2       834       315
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.027098e+02   8.88e-16
       4    9.625000e+01   2.62e-10   2.04e-07   2.70e-04      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   4.2% of solution time (0.000221 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.62500001106283e+01
Final feasibility error (abs / rel) =   2.62e-10 / 2.62e-10
Final optimality error  (abs / rel) =   2.04e-07 / 1.48e-08
# of iterations                     =          4 
# of CG iterations                  =          0 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00526 (     0.005 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  214 |                           214
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        76         0 |         2        76         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       460       329 |         2       460       329
  total:              2       840       329 |         2       840       329
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.030615e+02   1.78e-15
       5    9.800000e+01   8.86e-11   9.71e-10   2.67e-05      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   4.6% of solution time (0.000368 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.79999999985809e+01
Final feasibility error (abs / rel) =   8.86e-11 / 8.86e-11
Final optimality error  (abs / rel) =   9.71e-10 / 6.94e-11
# of iterations                     =          5 
# of CG iterations                  =          0 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00806 (     0.008 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  214 |                           214
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        76         0 |         2        76         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       460       331 |         2       460       331
  total:              2       840       331 |         2       840       331
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.044811e+02   8.88e-16
       7    9.308531e+01   6.60e-11   2.25e-08   1.18e-05      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   4.7% of solution time (0.000399 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.30853065879219e+01
Final feasibility error (abs / rel) =   6.60e-11 / 6.60e-11
Final optimality error  (abs / rel) =   2.25e-08 / 1.61e-09
# of iterations                     =          7 
# of CG iterations                  =          1 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00849 (     0.008 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  216 |                           216
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        78         0 |         2        78         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       472       333 |         2       472       333
  total:              2       852       333 |         2       852       333
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.045810e+02   1.78e-15
      10    9.148490e+01   1.79e-08   1.31e-05   1.69e-02      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   4.9% of solution time (0.000490 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.14848957551100e+01
Final feasibility error (abs / rel) =   1.79e-08 / 1.79e-08
Final optimality error  (abs / rel) =   1.31e-05 / 8.69e-07
# of iterations                     =         10 
# of CG iterations                  =          1 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.01001 (     0.010 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  210 |                           210
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        72         0 |         2        72         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       436       311 |         2       436       311
  total:              2       816       311 |         2       816       311
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.048688e+02   8.88e-16
       5    9.735387e+01   3.81e-07   6.83e-06   3.78e-03      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   6.2% of solution time (0.000482 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.73538727221633e+01
Final feasibility error (abs / rel) =   3.81e-07 / 3.81e-07
Final optimality error  (abs / rel) =   6.83e-06 / 4.91e-07
# of iterations                     =          5 
# of CG iterations                  =          3 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.00785 (     0.008 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  214 |                           214
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        76         0 |         2        76         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       460       321 |         2       460       321
  total:              2       840       321 |         2       840       321
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.049121e+02   8.88e-16
      10    9.953095e+01   4.11e-05   3.19e+00   5.80e-02      0.01
      20    9.943656e+01   1.61e-08   2.70e-02   2.58e-05      0.03
      23    9.943063e+01   5.09e-09   6.29e-06   4.20e-03      0.04

EXIT: Locally optimal solution found.

HINT: Knitro spent   1.1% of solution time (0.000456 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   9.94306346201769e+01
Final feasibility error (abs / rel) =   5.09e-09 / 5.09e-09
Final optimality error  (abs / rel) =   6.29e-06 / 5.06e-07
# of iterations                     =         23 
# of CG iterations                  =         28 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.03970 (     0.040 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  213 |                           213
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        75         0 |         2        75         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       454       323 |         2       454       323
  total:              2       834       323 |         2       834       323
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.059825e+02   1.78e-15
       8    8.830155e+01   8.70e-11   2.81e-09   3.16e-06      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   4.9% of solution time (0.000577 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   8.83015474514894e+01
Final feasibility error (abs / rel) =   8.70e-11 / 8.70e-11
Final optimality error  (abs / rel) =   2.81e-09 / 1.95e-10
# of iterations                     =          8 
# of CG iterations                  =          2 
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.01171 (     0.012 CPU time)
Time spent in evaluations (secs)    =       0.00000

================================================================================


=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro using 1 thread.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              0.001
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001

Problem Characteristics                     |           Presolved
-----------------------
Problem type: QCQP
Objective: minimize / quadratic
Number of variables:                     78 |                            78
  bounds:         lower     upper     range |     lower     upper     range
                     70         0         4 |        70         0         4
                             free     fixed |                free     fixed
                                4         0 |                   4         0
Number of constraints:                  217 |                           217
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            68        68         0 |        68        68         0
  quadratic:          2        79         0 |         2        79         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0       380           |         0       380          
  quadratic:          2       478       335 |         2       478       335
  total:              2       858       335 |         2       858       335
Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [1e+00, 8e+00] |                [1e+00, 8e+00]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 1e+00]

Knitro using the Interior-Point/Barrier Direct algorithm.

    Iter       Objective  FeasError   OptError   ||Step||      Time 
--------  --------------  ---------  ---------  ---------  --------
       0    1.071087e+02   1.78e-15
       3    1.040000e+02   1.79e-10   1.18e-05   2.44e-03      0.01

EXIT: Locally optimal solution found.

HINT: Knitro spent   8.0% of solution time (0.000474 secs) checking model
      convexity. To skip the automatic convexity checker for QPs and QCQPs,
      explicity set the user option convex=0 or convex=1.

Final Statistics
----------------
Final objective value               =   1.04000022642711e+02



          Commercial License
draw_areas(starts, packings)
Every start solvedstart 1: 65.85start 2: 89.71start 3: 91.55start 4: 94.63start 5: 95.67start 6: 96.33start 7: 98.22start 8: 98.78start 9: 100.99start 10: 101.16start 11: 101.51start 12: 101.56start 13: 102.71start 14: 103.06start 15: 104.48start 16: 104.58start 17: 104.87start 18: 104.91start 19: 105.98start 20: 107.11start 21: 107.11start 22: 107.77start 23: 107.93start 24: 108.24start 25: 109.40start 26: 110.10start 27: 110.36start 28: 110.63start 29: 112.13start 30: 112.94start 31: 113.58start 32: 113.86start 33: 118.31start 34: 119.51start 35: 120.34start 36: 121.37start 37: 122.62start 38: 129.28start 39: 130.01start 40: 132.26start 41: 133.79start 42: 133.94start 43: 134.20start 44: 134.63start 45: 135.64start 46: 136.42start 47: 137.25start 48: 137.49start 49: 137.71start 50: 137.76start 51: 137.94start 52: 138.09start 53: 139.13start 54: 139.22start 55: 140.01start 56: 141.76start 57: 144.15start 58: 144.44start 59: 144.78start 60: 145.08start 61: 145.38start 62: 145.64start 63: 148.90start 64: 153.51solved 1: 60.00solved 2: 88.65solved 3: 90.00solved 4: 90.00solved 5: 92.00solved 6: 93.31solved 7: 98.13solved 8: 93.19solved 9: 95.69solved 10: 100.47solved 11: 101.33solved 12: 95.23solved 13: 96.25solved 14: 98.00solved 15: 93.09solved 16: 91.48solved 17: 97.35solved 18: 99.43solved 19: 88.30solved 20: 104.00solved 21: 104.00solved 22: 93.98solved 23: 96.57solved 24: 97.61solved 25: 90.92solved 26: 104.13solved 27: 98.67solved 28: 97.06solved 29: 97.06solved 30: 100.27solved 31: 96.00solved 32: 97.33solved 33: 102.54solved 34: 117.70solved 35: 110.32solved 36: 106.67solved 37: 93.31solved 38: 95.40solved 39: 115.42solved 40: 96.24solved 41: 100.59solved 42: 112.53solved 43: 116.23solved 44: 100.23solved 45: 116.85solved 46: 101.81solved 47: 109.63solved 48: 104.71solved 49: 98.69solved 50: 96.13solved 51: 98.82solved 52: 115.74solved 53: 127.63solved 54: 99.01solved 55: 115.06solved 56: 96.13solved 57: 99.53solved 58: 97.23solved 59: 99.01solved 60: 113.46solved 61: 104.00solved 62: 101.16solved 63: 96.65solved 64: 104.00best 60.0010203040506050100150Start, ordered by increasing areaRectangle areastart areasolved area

Conclusion

report_spread(starts, packings)
     start   start area   solved area   of best
-----------------------------------------------
  smallest        65.85         60.00     1.00x
    median       118.31        102.54     1.71x
   largest       153.51        104.00     1.73x

64 starts solved, 1 reaching the best area of 60.00. The worst solve is 2.13 times larger.

The plot above shows the value of the initial point strategy: the model being nonconvex, the quality of the solutions degrades as the area of the initial points increases, and no initial point leads to a better solution than the best one.

Generating good initial points and solving the simplified model from them is therefore an effective strategy for this packing problem.