From a8794a8002bb07512e34cf967bbb094c0dd7366a Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Fri, 19 May 2023 12:00:18 -0600
Subject: [PATCH 01/15] farm_layout initial implementation

---
 examples/gdp/farm_layout/farm_layout.py | 118 ++++++++++++++++++++++++
 1 file changed, 118 insertions(+)
 create mode 100644 examples/gdp/farm_layout/farm_layout.py

diff --git a/examples/gdp/farm_layout/farm_layout.py b/examples/gdp/farm_layout/farm_layout.py
new file mode 100644
index 00000000000..ee984106bfd
--- /dev/null
+++ b/examples/gdp/farm_layout/farm_layout.py
@@ -0,0 +1,118 @@
+"""
+Farm layout example from Sawaya (2006)
+"""
+
+from __future__ import division
+
+from pyomo.environ import ConcreteModel, Objective, Param, RangeSet, Set, Var, Constraint, NonNegativeReals
+
+farm_layout_model_examples = {
+    "Flay02": [
+        {1: 40, 2: 50},
+        {1: 1, 2: 1},
+        {1: 1, 2: 1},
+        30,
+        30,
+    ]
+}
+
+def build_model(areas, length_lbs, width_lbs, length_upper_overall, width_upper_overall):
+    """Build the model."""
+    
+    # Ensure good data was passed
+    for a in areas:
+        assert a > 0, "Area must be positive"
+    
+    m = ConcreteModel(name="Farm Layout model")
+    m.plots = RangeSet(len(areas), doc=f"{len(areas)} plots")
+    m.plot_pairs = Set(
+        initialize=m.plots * m.plots, filter=lambda _, p1, p2: p1 < p2
+    )
+    
+    m.length_lbs = Param(m.plots, initialize=length_lbs, doc="Lower bounds for length of individual rectangles")
+    m.width_lbs = Param(m.plots, initialize=width_lbs, doc="Lower bounds for length of individual rectangles")
+    m.plot_area = Param(m.plots, initialize=areas, doc="Area of this plot")
+    
+    m.overall_length_ub = Param(initialize=length_upper_overall, doc="Overall upper bound for length")
+    m.overall_width_ub = Param(initialize=width_upper_overall, doc="Overall upper bound for width")
+
+    m.Length = Var(domain=NonNegativeReals, doc="Overall length we are minimizing")
+    m.Width = Var(domain=NonNegativeReals, doc="Overall width we are minimizing")
+    
+    m.perim = Objective(expr=2 * (m.Length + m.Width), doc="Perimeter of the farmland configuration")
+
+    m.plot_x = Var(m.plots, bounds=(0, m.overall_length_ub), doc="x-coordinate of plot center")
+    m.plot_y = Var(m.plots, bounds=(0, m.overall_width_ub), doc="y-coordinate of plot center")
+    m.plot_length = Var(m.plots, bounds=(0, m.overall_length_ub), doc="Length of this plot")
+    m.plot_width = Var(m.plots, bounds=(0, m.overall_width_ub), doc="Width of this plot")
+
+    # Constraints
+
+    # I think this constraint is actually wrong in the paper. Shouldn't
+    # it be plot_length / 2, since plot_x is the middle of the rectangle?
+    @m.Constraint(m.plots, doc="Length is consistent with lengths of plots")
+    def length_consistency(m, p):
+        return m.Length >= m.plot_x[p] + m.plot_length[p]
+    # As above 
+    @m.Constraint(m.plots, doc="Width is consistent with widths of plots")
+    def width_consistency(m, p):
+        return m.Width >= m.plot_y[p] + m.plot_width[p]
+    
+    @m.Constraint(m.plots, doc="Lengths and widths are consistent with area")
+    def area_consistency(m, p):
+        return m.plot_area[p] / m.plot_width[p] - m.plot_length[p] <= 0
+    
+    @m.Disjunction(
+        m.plot_pairs,
+        doc="Make sure that none of the plots overlap in "
+        "either the x or y dimensions.",
+    )
+    def no_overlap(m, p1, p2):
+        return [
+            m.plot_x[p1] + m.plot_length[p1] / 2
+            <= (m.plot_x[p2] - m.plot_length[p2] / 2),
+            m.plot_y[p1] + m.plot_width[p1] / 2
+            <= (m.plot_y[p2] - m.plot_width[p2] / 2),
+            m.plot_x[p2] + m.plot_length[p2] / 2
+            <= (m.plot_x[p1] - m.plot_length[p1] / 2),
+            m.plot_y[p2] + m.plot_width[p2] / 2
+            <= (m.plot_y[p1] - m.plot_width[p1] / 2),
+        ]
+    m.length_cons = Constraint(expr=m.Length <= m.overall_length_ub)
+    m.width_cons = Constraint(expr=m.Width <= m.overall_width_ub)
+
+    # what do these do exactly? Shouldn't this be dividing by width_lbs ??
+    @m.Constraint(m.plots, doc="Length bounds?")
+    def length_bounds(m, p):
+        return (m.length_lbs[p], m.plot_length[p], m.plot_area[p] / m.length_lbs[p])
+    @m.Constraint(m.plots, doc="Width bounds?")
+    def width_bounds(m, p):
+        return (m.width_lbs[p], m.plot_width[p],  m.plot_area[p] / m.width_lbs[p])
+
+    # ensure compatibility between coordinates, l/w lower bounds, and overall upper bounds.
+    # But I think this is also wrong. Shouldn't it use L_i^L / 2 ?
+    @m.Constraint(m.plots, doc="x-coordinate compatibility")
+    def x_compat(m, p):
+        return m.plot_x[p] <= m.overall_length_ub - m.length_lbs[p]
+    @m.Constraint(m.plots, doc="y-coordinate compatibility")
+    def y_compat(m, p):
+        return m.plot_y[p] <= m.overall_width_ub - m.width_lbs[p]
+
+    return m
+
+if __name__ == "__main__":
+    from pyomo.environ import SolverFactory
+    from pyomo.core.base import TransformationFactory
+
+    # Set up a solver, for example scip and bigm works
+    solver = SolverFactory("scip")
+    transformer = TransformationFactory("gdp.bigm")
+
+    print("building model")
+    model = build_model(*farm_layout_model_examples["Flay02"])
+    print("applying transformer")
+    transformer.apply_to(model)
+    print("solving example problem")
+    solver.solve(model)
+    print(f"Found objective function value: {model.perim()}")
+    print()

From d3abd49fbbeb192adf7b23e2a7adb3eb684e0add Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Fri, 19 May 2023 13:41:54 -0600
Subject: [PATCH 02/15] add drawing code

---
 examples/gdp/farm_layout/farm_layout.py | 56 +++++++++++++++++++++++++
 1 file changed, 56 insertions(+)

diff --git a/examples/gdp/farm_layout/farm_layout.py b/examples/gdp/farm_layout/farm_layout.py
index ee984106bfd..b010773cd65 100644
--- a/examples/gdp/farm_layout/farm_layout.py
+++ b/examples/gdp/farm_layout/farm_layout.py
@@ -100,6 +100,61 @@ def y_compat(m, p):
 
     return m
 
+
+def draw_model(m, title=None):
+    """Draw a model using matplotlib to illustrate what's going on. Pass 'title' kwarg to give chart a title"""
+
+    # matplotlib setup
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+
+    fig, ax = plt.subplots(1, 1, figsize=(10, 10))
+    ax.set_xlabel("x")
+    ax.set_ylabel("y")
+
+    # Set bounds
+    ax.set_xlim([m.overall_length_ub() * -0.2, m.overall_length_ub() * 1.2 ])
+    ax.set_ylim([m.overall_width_ub() * -0.2, m.overall_width_ub() * 1.2])
+
+    if title is not None:
+        plt.title(title)
+    
+    # First, let's transparently draw the overall bounds
+    ax.add_patch(
+        mpl.patches.Rectangle(
+        (0, 0),
+        m.overall_length_ub(),
+        m.overall_width_ub(),
+        facecolor="#229922",
+        edgecolor="#000000",
+        alpha=0.2
+        )
+    )
+
+    for p in m.plots:
+        print(f"drawing plot {p}: x={m.plot_x[p]()}, y={m.plot_y[p]()}, length={m.plot_length[p]()}, width={m.plot_width[p]()}")
+        ax.add_patch(
+            mpl.patches.Rectangle(
+            (
+                m.plot_x[p]() - m.plot_length[p]() / 2,
+                m.plot_y[p]() - m.plot_width[p]() / 2,
+            ),
+            m.plot_length[p](),
+            m.plot_width[p](),
+            facecolor="#ebdc78",
+            edgecolor="#000000",
+            )
+        )
+        ax.text(
+            m.plot_x[p](),
+            m.plot_y[p](),
+            f"Plot {p}",
+            horizontalalignment="center",
+            verticalalignment="center",
+        )
+
+    plt.show()
+
 if __name__ == "__main__":
     from pyomo.environ import SolverFactory
     from pyomo.core.base import TransformationFactory
@@ -115,4 +170,5 @@ def y_compat(m, p):
     print("solving example problem")
     solver.solve(model)
     print(f"Found objective function value: {model.perim()}")
+    draw_model(model)
     print()

From 359416d63296bb8762283f62c197996e9e6df7ae Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Tue, 23 May 2023 14:51:29 -0600
Subject: [PATCH 03/15] fix model for sane-ness in contradiction to paper

Now objective values match those generated by scip using the paper's published .nl files. Also add the rest of the instances.
---
 examples/gdp/farm_layout/farm_layout.py | 117 ++++++++++++++++--------
 1 file changed, 79 insertions(+), 38 deletions(-)

diff --git a/examples/gdp/farm_layout/farm_layout.py b/examples/gdp/farm_layout/farm_layout.py
index b010773cd65..1d6324989ce 100644
--- a/examples/gdp/farm_layout/farm_layout.py
+++ b/examples/gdp/farm_layout/farm_layout.py
@@ -1,25 +1,61 @@
 """
-Farm layout example from Sawaya (2006)
+Farm layout example from Sawaya (2006). This is a GDP problem with some hyperbolic constraints to establish
+consistency of areas with length and width. The FLay05 and FLay06 instances may take some time to solve; the
+others should be fast. Note that the Sawaya paper contains a little bit of nonclarity: it references "height"
+variables which do not exist - we use "length" for the x-axis and "width" on the y-axis, and it also is
+unclear on the way the coordinates define the rectangles; we have decided that they are on the bottom-left
+and adapted the disjunction constraints to match.
 """
 
 from __future__ import division
 
 from pyomo.environ import ConcreteModel, Objective, Param, RangeSet, Set, Var, Constraint, NonNegativeReals
 
+# Format: areas, length lower bounds, width lower bounds, length overall upper bound, width overall upper bound
 farm_layout_model_examples = {
-    "Flay02": [
+    "FLay02": [
         {1: 40, 2: 50},
         {1: 1, 2: 1},
         {1: 1, 2: 1},
         30,
         30,
+    ],
+    "FLay03": [
+        {1: 40, 2: 50, 3: 60},
+        {1: 1, 2: 1, 3: 1},
+        {1: 1, 2: 1, 3: 1},
+        30,
+        30,
+    ],
+    "FLay04": [
+        {1: 40, 2: 50, 3: 60, 4: 35},
+        {1: 3, 2: 3, 3: 3, 4: 3},
+        {1: 3, 2: 3, 3: 3, 4: 3},
+        100,
+        100,
+    ],
+    "FLay05": [
+        {1: 40, 2: 50, 3: 60, 4: 35, 5: 75},
+        {1: 1, 2: 1, 3: 1, 4: 1, 5: 1},
+        {1: 1, 2: 1, 3: 1, 4: 1, 5: 1},
+        30,
+        30,
+    ],
+    "FLay06": [
+        {1: 40, 2: 50, 3: 60, 4: 35, 5: 75, 6: 20},
+        {1: 1, 2: 1, 3: 1, 4: 1, 5: 1, 6: 1},
+        {1: 1, 2: 1, 3: 1, 4: 1, 5: 1, 6: 1},
+        30,
+        30,
     ]
 }
 
-def build_model(areas, length_lbs, width_lbs, length_upper_overall, width_upper_overall):
-    """Build the model."""
+def build_model(params=farm_layout_model_examples["FLay03"]):
+    """Build the model. Default example gives FLay03 which is fairly small"""
     
     # Ensure good data was passed
+    assert len(params) == 5, "Params should be in the format: areas, length_lbs, width_lbs, length_upper_overall, width_upper, overall"
+    areas, length_lbs, width_lbs, length_upper_overall, width_upper_overall = params
     for a in areas:
         assert a > 0, "Area must be positive"
     
@@ -36,8 +72,8 @@ def build_model(areas, length_lbs, width_lbs, length_upper_overall, width_upper_
     m.overall_length_ub = Param(initialize=length_upper_overall, doc="Overall upper bound for length")
     m.overall_width_ub = Param(initialize=width_upper_overall, doc="Overall upper bound for width")
 
-    m.Length = Var(domain=NonNegativeReals, doc="Overall length we are minimizing")
-    m.Width = Var(domain=NonNegativeReals, doc="Overall width we are minimizing")
+    m.Length = Var(domain=NonNegativeReals, doc="Overall length of the farmland configuration")
+    m.Width = Var(domain=NonNegativeReals, doc="Overall width of the farmland configuration")
     
     m.perim = Objective(expr=2 * (m.Length + m.Width), doc="Perimeter of the farmland configuration")
 
@@ -47,9 +83,8 @@ def build_model(areas, length_lbs, width_lbs, length_upper_overall, width_upper_
     m.plot_width = Var(m.plots, bounds=(0, m.overall_width_ub), doc="Width of this plot")
 
     # Constraints
-
-    # I think this constraint is actually wrong in the paper. Shouldn't
-    # it be plot_length / 2, since plot_x is the middle of the rectangle?
+    
+    # Ensure consistency of Length with lengths and x-values of plots
     @m.Constraint(m.plots, doc="Length is consistent with lengths of plots")
     def length_consistency(m, p):
         return m.Length >= m.plot_x[p] + m.plot_length[p]
@@ -69,28 +104,28 @@ def area_consistency(m, p):
     )
     def no_overlap(m, p1, p2):
         return [
-            m.plot_x[p1] + m.plot_length[p1] / 2
-            <= (m.plot_x[p2] - m.plot_length[p2] / 2),
-            m.plot_y[p1] + m.plot_width[p1] / 2
-            <= (m.plot_y[p2] - m.plot_width[p2] / 2),
-            m.plot_x[p2] + m.plot_length[p2] / 2
-            <= (m.plot_x[p1] - m.plot_length[p1] / 2),
-            m.plot_y[p2] + m.plot_width[p2] / 2
-            <= (m.plot_y[p1] - m.plot_width[p1] / 2),
+            m.plot_x[p1] + m.plot_length[p1]
+            <= m.plot_x[p2],
+            m.plot_y[p1] + m.plot_width[p1]
+            <= m.plot_y[p2],
+            m.plot_x[p2] + m.plot_length[p2]
+            <= m.plot_x[p1],
+            m.plot_y[p2] + m.plot_width[p2]
+            <= m.plot_y[p1] 
         ]
     m.length_cons = Constraint(expr=m.Length <= m.overall_length_ub)
     m.width_cons = Constraint(expr=m.Width <= m.overall_width_ub)
 
-    # what do these do exactly? Shouldn't this be dividing by width_lbs ??
-    @m.Constraint(m.plots, doc="Length bounds?")
+    # This imposes a square-ness constraint. I'm not sure the justification for these 
+    # upper bounds in particular
+    @m.Constraint(m.plots, doc="Length bounds")
     def length_bounds(m, p):
         return (m.length_lbs[p], m.plot_length[p], m.plot_area[p] / m.length_lbs[p])
-    @m.Constraint(m.plots, doc="Width bounds?")
+    @m.Constraint(m.plots, doc="Width bounds")
     def width_bounds(m, p):
         return (m.width_lbs[p], m.plot_width[p],  m.plot_area[p] / m.width_lbs[p])
 
     # ensure compatibility between coordinates, l/w lower bounds, and overall upper bounds.
-    # But I think this is also wrong. Shouldn't it use L_i^L / 2 ?
     @m.Constraint(m.plots, doc="x-coordinate compatibility")
     def x_compat(m, p):
         return m.plot_x[p] <= m.overall_length_ub - m.length_lbs[p]
@@ -102,7 +137,7 @@ def y_compat(m, p):
 
 
 def draw_model(m, title=None):
-    """Draw a model using matplotlib to illustrate what's going on. Pass 'title' kwarg to give chart a title"""
+    """Draw a model using matplotlib to illustrate what's going on. Pass 'title' arg to give chart a title"""
 
     # matplotlib setup
     import matplotlib as mpl
@@ -131,13 +166,14 @@ def draw_model(m, title=None):
         )
     )
 
+    # Now draw the plots
     for p in m.plots:
         print(f"drawing plot {p}: x={m.plot_x[p]()}, y={m.plot_y[p]()}, length={m.plot_length[p]()}, width={m.plot_width[p]()}")
         ax.add_patch(
             mpl.patches.Rectangle(
             (
-                m.plot_x[p]() - m.plot_length[p]() / 2,
-                m.plot_y[p]() - m.plot_width[p]() / 2,
+                m.plot_x[p](),
+                m.plot_y[p]()
             ),
             m.plot_length[p](),
             m.plot_width[p](),
@@ -146,8 +182,8 @@ def draw_model(m, title=None):
             )
         )
         ax.text(
-            m.plot_x[p](),
-            m.plot_y[p](),
+            m.plot_x[p]() + m.plot_length[p]() / 2,
+            m.plot_y[p]() + m.plot_width[p]() / 2,
             f"Plot {p}",
             horizontalalignment="center",
             verticalalignment="center",
@@ -159,16 +195,21 @@ def draw_model(m, title=None):
     from pyomo.environ import SolverFactory
     from pyomo.core.base import TransformationFactory
 
-    # Set up a solver, for example scip and bigm works
+    # Set up a solver, for example scip and hull works.
+    # Note that these constraints are not linear or quadratic.
     solver = SolverFactory("scip")
-    transformer = TransformationFactory("gdp.bigm")
-
-    print("building model")
-    model = build_model(*farm_layout_model_examples["Flay02"])
-    print("applying transformer")
-    transformer.apply_to(model)
-    print("solving example problem")
-    solver.solve(model)
-    print(f"Found objective function value: {model.perim()}")
-    draw_model(model)
-    print()
+    transformer = TransformationFactory("gdp.hull")
+
+    # Try all the instances except FLay05 and FLay06, since they take a while.
+    for key in farm_layout_model_examples.keys():
+        if key not in ["FLay06", "FLay05"]:
+            print(f"solving example problem: {key}")
+            print("building model")
+            model = build_model(farm_layout_model_examples[key])
+            print("applying transformer")
+            transformer.apply_to(model)
+            print("solving model")
+            solver.solve(model)
+            print(f"Found objective function value: {model.perim()}")
+            draw_model(model, title=key)
+            print()

From b7ee41521ed63ec5e7e150d3f99d7b4e26996acb Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Tue, 23 May 2023 14:57:15 -0600
Subject: [PATCH 04/15] black reformatting

---
 examples/gdp/farm_layout/farm_layout.py | 170 ++++++++++++++----------
 1 file changed, 98 insertions(+), 72 deletions(-)

diff --git a/examples/gdp/farm_layout/farm_layout.py b/examples/gdp/farm_layout/farm_layout.py
index 1d6324989ce..6aa9e10fc58 100644
--- a/examples/gdp/farm_layout/farm_layout.py
+++ b/examples/gdp/farm_layout/farm_layout.py
@@ -1,32 +1,31 @@
 """
-Farm layout example from Sawaya (2006). This is a GDP problem with some hyperbolic constraints to establish
-consistency of areas with length and width. The FLay05 and FLay06 instances may take some time to solve; the
-others should be fast. Note that the Sawaya paper contains a little bit of nonclarity: it references "height"
-variables which do not exist - we use "length" for the x-axis and "width" on the y-axis, and it also is
-unclear on the way the coordinates define the rectangles; we have decided that they are on the bottom-left
-and adapted the disjunction constraints to match.
+Farm layout example from Sawaya (2006). The goal is to determine optimal placements and dimensions for farm
+plots of specified areas to minimize the perimeter of a minimal enclosing fence. This is a GDP problem with 
+some hyperbolic constraints to establish consistency of areas with length and width. The FLay05 and FLay06 
+instances may take some time to solve; the others should be fast. Note that the Sawaya paper contains a 
+little bit of nonclarity: it references "height" variables which do not exist - we use "length" for the x-axis 
+and "width" on the y-axis, and it also is unclear on the way the coordinates define the rectangles; we have
+ decided that they are on the bottom-left and adapted the disjunction constraints to match.
 """
 
 from __future__ import division
 
-from pyomo.environ import ConcreteModel, Objective, Param, RangeSet, Set, Var, Constraint, NonNegativeReals
+from pyomo.environ import (
+    ConcreteModel,
+    Objective,
+    Param,
+    RangeSet,
+    Set,
+    Var,
+    Constraint,
+    NonNegativeReals,
+)
 
 # Format: areas, length lower bounds, width lower bounds, length overall upper bound, width overall upper bound
+# Examples are from Sawaya (2006)
 farm_layout_model_examples = {
-    "FLay02": [
-        {1: 40, 2: 50},
-        {1: 1, 2: 1},
-        {1: 1, 2: 1},
-        30,
-        30,
-    ],
-    "FLay03": [
-        {1: 40, 2: 50, 3: 60},
-        {1: 1, 2: 1, 3: 1},
-        {1: 1, 2: 1, 3: 1},
-        30,
-        30,
-    ],
+    "FLay02": [{1: 40, 2: 50}, {1: 1, 2: 1}, {1: 1, 2: 1}, 30, 30],
+    "FLay03": [{1: 40, 2: 50, 3: 60}, {1: 1, 2: 1, 3: 1}, {1: 1, 2: 1, 3: 1}, 30, 30],
     "FLay04": [
         {1: 40, 2: 50, 3: 60, 4: 35},
         {1: 3, 2: 3, 3: 3, 4: 3},
@@ -47,56 +46,84 @@
         {1: 1, 2: 1, 3: 1, 4: 1, 5: 1, 6: 1},
         30,
         30,
-    ]
+    ],
 }
 
+
 def build_model(params=farm_layout_model_examples["FLay03"]):
     """Build the model. Default example gives FLay03 which is fairly small"""
-    
+
     # Ensure good data was passed
-    assert len(params) == 5, "Params should be in the format: areas, length_lbs, width_lbs, length_upper_overall, width_upper, overall"
+    assert (
+        len(params) == 5
+    ), "Params should be in the format: areas, length_lbs, width_lbs, length_upper_overall, width_upper, overall"
     areas, length_lbs, width_lbs, length_upper_overall, width_upper_overall = params
     for a in areas:
         assert a > 0, "Area must be positive"
-    
+
     m = ConcreteModel(name="Farm Layout model")
     m.plots = RangeSet(len(areas), doc=f"{len(areas)} plots")
-    m.plot_pairs = Set(
-        initialize=m.plots * m.plots, filter=lambda _, p1, p2: p1 < p2
+    m.plot_pairs = Set(initialize=m.plots * m.plots, filter=lambda _, p1, p2: p1 < p2)
+
+    m.length_lbs = Param(
+        m.plots,
+        initialize=length_lbs,
+        doc="Lower bounds for length of individual rectangles",
+    )
+    m.width_lbs = Param(
+        m.plots,
+        initialize=width_lbs,
+        doc="Lower bounds for length of individual rectangles",
     )
-    
-    m.length_lbs = Param(m.plots, initialize=length_lbs, doc="Lower bounds for length of individual rectangles")
-    m.width_lbs = Param(m.plots, initialize=width_lbs, doc="Lower bounds for length of individual rectangles")
     m.plot_area = Param(m.plots, initialize=areas, doc="Area of this plot")
-    
-    m.overall_length_ub = Param(initialize=length_upper_overall, doc="Overall upper bound for length")
-    m.overall_width_ub = Param(initialize=width_upper_overall, doc="Overall upper bound for width")
 
-    m.Length = Var(domain=NonNegativeReals, doc="Overall length of the farmland configuration")
-    m.Width = Var(domain=NonNegativeReals, doc="Overall width of the farmland configuration")
-    
-    m.perim = Objective(expr=2 * (m.Length + m.Width), doc="Perimeter of the farmland configuration")
+    m.overall_length_ub = Param(
+        initialize=length_upper_overall, doc="Overall upper bound for length"
+    )
+    m.overall_width_ub = Param(
+        initialize=width_upper_overall, doc="Overall upper bound for width"
+    )
 
-    m.plot_x = Var(m.plots, bounds=(0, m.overall_length_ub), doc="x-coordinate of plot center")
-    m.plot_y = Var(m.plots, bounds=(0, m.overall_width_ub), doc="y-coordinate of plot center")
-    m.plot_length = Var(m.plots, bounds=(0, m.overall_length_ub), doc="Length of this plot")
-    m.plot_width = Var(m.plots, bounds=(0, m.overall_width_ub), doc="Width of this plot")
+    m.Length = Var(
+        domain=NonNegativeReals, doc="Overall length of the farmland configuration"
+    )
+    m.Width = Var(
+        domain=NonNegativeReals, doc="Overall width of the farmland configuration"
+    )
+
+    m.perim = Objective(
+        expr=2 * (m.Length + m.Width), doc="Perimeter of the farmland configuration"
+    )
+
+    m.plot_x = Var(
+        m.plots, bounds=(0, m.overall_length_ub), doc="x-coordinate of plot center"
+    )
+    m.plot_y = Var(
+        m.plots, bounds=(0, m.overall_width_ub), doc="y-coordinate of plot center"
+    )
+    m.plot_length = Var(
+        m.plots, bounds=(0, m.overall_length_ub), doc="Length of this plot"
+    )
+    m.plot_width = Var(
+        m.plots, bounds=(0, m.overall_width_ub), doc="Width of this plot"
+    )
 
     # Constraints
-    
+
     # Ensure consistency of Length with lengths and x-values of plots
     @m.Constraint(m.plots, doc="Length is consistent with lengths of plots")
     def length_consistency(m, p):
         return m.Length >= m.plot_x[p] + m.plot_length[p]
-    # As above 
+
+    # As above
     @m.Constraint(m.plots, doc="Width is consistent with widths of plots")
     def width_consistency(m, p):
         return m.Width >= m.plot_y[p] + m.plot_width[p]
-    
+
     @m.Constraint(m.plots, doc="Lengths and widths are consistent with area")
     def area_consistency(m, p):
         return m.plot_area[p] / m.plot_width[p] - m.plot_length[p] <= 0
-    
+
     @m.Disjunction(
         m.plot_pairs,
         doc="Make sure that none of the plots overlap in "
@@ -104,31 +131,30 @@ def area_consistency(m, p):
     )
     def no_overlap(m, p1, p2):
         return [
-            m.plot_x[p1] + m.plot_length[p1]
-            <= m.plot_x[p2],
-            m.plot_y[p1] + m.plot_width[p1]
-            <= m.plot_y[p2],
-            m.plot_x[p2] + m.plot_length[p2]
-            <= m.plot_x[p1],
-            m.plot_y[p2] + m.plot_width[p2]
-            <= m.plot_y[p1] 
+            m.plot_x[p1] + m.plot_length[p1] <= m.plot_x[p2],
+            m.plot_y[p1] + m.plot_width[p1] <= m.plot_y[p2],
+            m.plot_x[p2] + m.plot_length[p2] <= m.plot_x[p1],
+            m.plot_y[p2] + m.plot_width[p2] <= m.plot_y[p1],
         ]
+
     m.length_cons = Constraint(expr=m.Length <= m.overall_length_ub)
     m.width_cons = Constraint(expr=m.Width <= m.overall_width_ub)
 
-    # This imposes a square-ness constraint. I'm not sure the justification for these 
+    # This imposes a square-ness constraint. I'm not sure the justification for these
     # upper bounds in particular
     @m.Constraint(m.plots, doc="Length bounds")
     def length_bounds(m, p):
         return (m.length_lbs[p], m.plot_length[p], m.plot_area[p] / m.length_lbs[p])
+
     @m.Constraint(m.plots, doc="Width bounds")
     def width_bounds(m, p):
-        return (m.width_lbs[p], m.plot_width[p],  m.plot_area[p] / m.width_lbs[p])
+        return (m.width_lbs[p], m.plot_width[p], m.plot_area[p] / m.width_lbs[p])
 
     # ensure compatibility between coordinates, l/w lower bounds, and overall upper bounds.
     @m.Constraint(m.plots, doc="x-coordinate compatibility")
     def x_compat(m, p):
         return m.plot_x[p] <= m.overall_length_ub - m.length_lbs[p]
+
     @m.Constraint(m.plots, doc="y-coordinate compatibility")
     def y_compat(m, p):
         return m.plot_y[p] <= m.overall_width_ub - m.width_lbs[p]
@@ -148,37 +174,36 @@ def draw_model(m, title=None):
     ax.set_ylabel("y")
 
     # Set bounds
-    ax.set_xlim([m.overall_length_ub() * -0.2, m.overall_length_ub() * 1.2 ])
+    ax.set_xlim([m.overall_length_ub() * -0.2, m.overall_length_ub() * 1.2])
     ax.set_ylim([m.overall_width_ub() * -0.2, m.overall_width_ub() * 1.2])
 
     if title is not None:
         plt.title(title)
-    
+
     # First, let's transparently draw the overall bounds
     ax.add_patch(
         mpl.patches.Rectangle(
-        (0, 0),
-        m.overall_length_ub(),
-        m.overall_width_ub(),
-        facecolor="#229922",
-        edgecolor="#000000",
-        alpha=0.2
+            (0, 0),
+            m.overall_length_ub(),
+            m.overall_width_ub(),
+            facecolor="#229922",
+            edgecolor="#000000",
+            alpha=0.2,
         )
     )
 
     # Now draw the plots
     for p in m.plots:
-        print(f"drawing plot {p}: x={m.plot_x[p]()}, y={m.plot_y[p]()}, length={m.plot_length[p]()}, width={m.plot_width[p]()}")
+        print(
+            f"drawing plot {p}: x={m.plot_x[p]()}, y={m.plot_y[p]()}, length={m.plot_length[p]()}, width={m.plot_width[p]()}"
+        )
         ax.add_patch(
             mpl.patches.Rectangle(
-            (
-                m.plot_x[p](),
-                m.plot_y[p]()
-            ),
-            m.plot_length[p](),
-            m.plot_width[p](),
-            facecolor="#ebdc78",
-            edgecolor="#000000",
+                (m.plot_x[p](), m.plot_y[p]()),
+                m.plot_length[p](),
+                m.plot_width[p](),
+                facecolor="#ebdc78",
+                edgecolor="#000000",
             )
         )
         ax.text(
@@ -191,6 +216,7 @@ def draw_model(m, title=None):
 
     plt.show()
 
+
 if __name__ == "__main__":
     from pyomo.environ import SolverFactory
     from pyomo.core.base import TransformationFactory

From 504aa9bf8f7c750df5410cdb2990d13313e53972 Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Tue, 23 May 2023 14:58:02 -0600
Subject: [PATCH 05/15] fix a character

---
 examples/gdp/farm_layout/farm_layout.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/gdp/farm_layout/farm_layout.py b/examples/gdp/farm_layout/farm_layout.py
index 6aa9e10fc58..f85beaf5af1 100644
--- a/examples/gdp/farm_layout/farm_layout.py
+++ b/examples/gdp/farm_layout/farm_layout.py
@@ -5,7 +5,7 @@
 instances may take some time to solve; the others should be fast. Note that the Sawaya paper contains a 
 little bit of nonclarity: it references "height" variables which do not exist - we use "length" for the x-axis 
 and "width" on the y-axis, and it also is unclear on the way the coordinates define the rectangles; we have
- decided that they are on the bottom-left and adapted the disjunction constraints to match.
+decided that they are on the bottom-left and adapted the disjunction constraints to match.
 """
 
 from __future__ import division

From 2229923cca51369e54ca64afa87c408019b86458 Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Tue, 30 May 2023 13:50:42 -0600
Subject: [PATCH 06/15] add some alternate examples that prevent the solution
 from simply being a square

---
 examples/gdp/farm_layout/farm_layout.py | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/examples/gdp/farm_layout/farm_layout.py b/examples/gdp/farm_layout/farm_layout.py
index f85beaf5af1..0fc071df90a 100644
--- a/examples/gdp/farm_layout/farm_layout.py
+++ b/examples/gdp/farm_layout/farm_layout.py
@@ -47,6 +47,21 @@
         30,
         30,
     ],
+    # An example where the dimension constraints are strict enough that it cannot fit a square
+    "FLay03_alt_1": [
+        {1: 40, 2: 50, 3: 60},
+        {1: 5, 2: 5, 3: 5},
+        {1: 5, 2: 5, 3: 7},
+        30,
+        30,
+    ],
+    "FLay03_alt_2": [
+        {1: 40, 2: 50, 3: 60},
+        {1: 1, 2: 1, 3: 1},
+        {1: 1, 2: 1, 3: 1},
+        8,
+        30,
+    ],
 }
 
 

From 39549a55adc8378007e88634549f47838e61bc19 Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Tue, 30 May 2023 12:15:15 -0600
Subject: [PATCH 07/15] Add circles gdp example

---
 examples/gdp/circles/circles.py | 201 ++++++++++++++++++++++++++++++++
 1 file changed, 201 insertions(+)
 create mode 100644 examples/gdp/circles/circles.py

diff --git a/examples/gdp/circles/circles.py b/examples/gdp/circles/circles.py
new file mode 100644
index 00000000000..48ae13a280f
--- /dev/null
+++ b/examples/gdp/circles/circles.py
@@ -0,0 +1,201 @@
+"""
+The "circles" GDP example problem originating in Lee and Grossman (2000). The 
+goal is to choose a point to minimize a convex quadratic function over a set of 
+disjoint hyperspheres.
+"""
+
+from __future__ import division
+
+from pyomo.environ import (
+    ConcreteModel,
+    Objective,
+    Param,
+    RangeSet,
+    Var,
+    Constraint,
+    value,
+)
+
+# Circles2D3 is the original example described in the papers. Ruiz and Grossman 
+# (2012) discuss two more example instances, Circles2D36 and Circles3D36, but I 
+# couldn't find them.
+# format: dimension, circle_centers, circle_rvals, circle_penalties, reference_point, upper_bound
+circles_model_examples = {
+    "Circles2D3": [
+        2,
+        #{1: [0, 0], 2: [4, 1], 3: [2, 4]},
+        {(1,1): 0, (1,2): 0,
+         (2,1): 4, (2,2): 1,
+         (3,1): 2, (3,2): 4},
+        {1: 1, 2: 1, 3: 1},
+        {1: 2, 2: 1, 3: 3},
+        {1: 3, 2: 2},
+        8
+    ],
+    # Here the solver has a local minimum to avoid by not choosing the closest point
+    "Circles2D3_modified": [
+        2,
+        {(1,1): 0, (1,2): 0,
+         (2,1): 4, (2,2): 1,
+         (3,1): 2, (3,2): 4},
+        {1: 1, 2: 1, 3: 1},
+        {1: 2, 2: 3, 3: 1},
+        {1: 3, 2: 2},
+        8
+    ],
+    # Here's a 3D model.
+    "Circles3D4": [
+        3,
+        {
+            (1, 1): 0, (1, 2): 0, (1, 3): 0,
+            (2, 1): 2, (2, 2): 0, (2, 3): 1,
+            (3, 1): 3, (3, 2): 3, (3, 3): 3,
+            (4, 1): 1, (4, 2): 6, (4, 3): 1,
+        },
+        {1: 1, 2: 1, 3: 1, 4: 1},
+        {1: 1, 2: 3, 3: 2, 4: 1},
+        {1: 2, 2: 2, 3: 1},
+        10
+    ]
+}
+
+def build_model(data):
+    """Build a model. By default you get Circles2D3, a small instance."""
+
+    # Ensure good data was passed
+    assert len(data) == 6, "Error processing data"
+    dimension, circle_centers, circle_rvals, circle_penalties, reference_point, upper_bound = data
+
+    assert len(circle_rvals) == len(circle_penalties), "Error processing data"
+    assert len(circle_centers) == len(circle_rvals) * dimension, "Error processing data"
+    assert len(reference_point) == dimension, "Error processing data"
+
+    m = ConcreteModel()
+
+    m.circles = RangeSet(len(circle_rvals), doc=f"{len(circle_rvals)} circles")
+    m.idx = RangeSet(dimension, doc="n-dimensional indexing set for coordinates")
+
+    m.circ_centers = Param(m.circles, m.idx, initialize=circle_centers, doc="Center points of the circles")
+    m.circ_rvals = Param(m.circles, initialize=circle_rvals, doc="Squared radii of circles")
+    m.ref_point = Param(m.idx, initialize=reference_point, doc="Reference point for distance calculations")
+    m.circ_penalties = Param(m.circles, initialize=circle_penalties, doc="Penalty for being in each circle")
+
+    # Choose a point to minimize objective
+    m.point = Var(m.idx, bounds=(0, upper_bound), doc="Chosen point at which we evaluate objective function")
+
+    # Let's set up the "active penalty" this way
+    m.active_penalty = Var(bounds=(0, max(m.circ_penalties[i] for i in m.circles)), doc="Penalty for being in the current circle")
+
+    # Disjunction: we must be in at least one circle (in fact exactly one since they are disjoint)
+    @m.Disjunct(m.circles)
+    def circle_disjunct(d, circ):
+        m = d.model()
+        d.inside_circle = Constraint(expr=sum((m.point[i] - m.circ_centers[circ, i]) ** 2 for i in m.idx) <= m.circ_rvals[circ])
+        d.penalty = Constraint(expr=m.active_penalty == m.circ_penalties[circ])
+    @m.Disjunction()
+    def circle_disjunction(m):
+        return [m.circle_disjunct[i] for i in m.circles]
+
+    # Objective function is Euclidean distance from reference point, plus a penalty depending on which circle we are in
+    m.cost = Objective(expr=sum((m.point[i] - m.ref_point[i]) ** 2 for i in m.idx) + m.active_penalty, doc="Distance from reference point plus penalty")   
+
+    return m
+
+def draw_model(m, title=None):
+    """Draw a model using matplotlib to illustrate what's going on. Pass 'title' arg to give chart a title"""
+
+    if len(m.idx) != 2:
+        print(f"Unable to draw: this drawing code only supports 2D models, but a {len(m.idx)}-dimensional model was passed.")
+        return
+
+    # matplotlib setup
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+
+    fig, ax = plt.subplots(1, 1, figsize=(10, 10))
+    ax.set_xlabel("x")
+    ax.set_ylabel("y")
+
+    ax.set_xlim([0, 10])
+    ax.set_ylim([0, 10])
+
+    if title is not None:
+        plt.title(title)
+
+    for c in m.circles:
+        x = m.circ_centers[c, 1]
+        y = m.circ_centers[c, 2]
+        # remember, we are keeping around the *squared* radii for some reason
+        r = m.circ_rvals[c] ** 0.5
+        penalty = m.circ_penalties[c]
+        print(f"drawing circle {c}: x={x}, y={y}, r={r}, penalty={penalty}")
+
+        ax.add_patch(
+            mpl.patches.Circle(
+                (x, y),
+                radius=r,
+                facecolor="#0d98e6",
+                edgecolor="#000000",
+            )
+        )
+        # the alignment appears to work by magic
+        ax.text(
+            x,
+            y + r + 0.3,
+            f"Circle {c}: penalty {penalty}",
+            horizontalalignment="center",
+        )
+
+    # now the reference point
+    ax.add_patch(
+        mpl.patches.Circle(
+            (m.ref_point[1], m.ref_point[2]),
+            radius=0.05,
+            facecolor="#902222",
+            edgecolor="#902222",
+        )
+    )
+    ax.text(
+        m.ref_point[1],
+        m.ref_point[2] + 0.15,
+        "Reference Point",
+        horizontalalignment="center",
+    )
+
+    # and the chosen point
+    ax.add_patch(
+        mpl.patches.Circle(
+            (value(m.point[1]), value(m.point[2])),
+            radius=0.05,
+            facecolor="#11BB11",
+            edgecolor="#11BB11",
+        )
+    )
+    ax.text(
+        value(m.point[1]),
+        value(m.point[2]) + 0.15,
+        "Chosen Point",
+        horizontalalignment="center",
+    )
+
+    plt.show()
+
+
+
+if __name__ == "__main__":
+    from pyomo.environ import SolverFactory
+    from pyomo.core.base import TransformationFactory
+
+    # Set up a solver, for example scip and bigm works
+    solver = SolverFactory("scip")
+    transformer = TransformationFactory("gdp.bigm")
+    
+
+    for key in circles_model_examples:
+        model = build_model(circles_model_examples[key])
+        print(f"Solving model {key}")
+        transformer.apply_to(model)
+        solver.solve(model)
+        pt_string ="(" + "".join(str(value(model.point[i])) + (", " if i != len(model.idx) else "") for i in model.idx) + ")"
+        print(f"Optimal value found: {model.cost()} with chosen point {pt_string}")
+        draw_model(model, title=key)
\ No newline at end of file

From 523b4e44e8e808b22970bd15d4d1e24f738c40dc Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Tue, 30 May 2023 13:53:10 -0600
Subject: [PATCH 08/15] apply black

---
 examples/gdp/circles/circles.py | 112 ++++++++++++++++++++++----------
 1 file changed, 77 insertions(+), 35 deletions(-)

diff --git a/examples/gdp/circles/circles.py b/examples/gdp/circles/circles.py
index 48ae13a280f..74cc74e7957 100644
--- a/examples/gdp/circles/circles.py
+++ b/examples/gdp/circles/circles.py
@@ -16,55 +16,67 @@
     value,
 )
 
-# Circles2D3 is the original example described in the papers. Ruiz and Grossman 
-# (2012) discuss two more example instances, Circles2D36 and Circles3D36, but I 
+# Circles2D3 is the original example described in the papers. Ruiz and Grossman
+# (2012) discuss two more example instances, Circles2D36 and Circles3D36, but I
 # couldn't find them.
 # format: dimension, circle_centers, circle_rvals, circle_penalties, reference_point, upper_bound
 circles_model_examples = {
     "Circles2D3": [
         2,
-        #{1: [0, 0], 2: [4, 1], 3: [2, 4]},
-        {(1,1): 0, (1,2): 0,
-         (2,1): 4, (2,2): 1,
-         (3,1): 2, (3,2): 4},
+        # {1: [0, 0], 2: [4, 1], 3: [2, 4]},
+        {(1, 1): 0, (1, 2): 0, (2, 1): 4, (2, 2): 1, (3, 1): 2, (3, 2): 4},
         {1: 1, 2: 1, 3: 1},
         {1: 2, 2: 1, 3: 3},
         {1: 3, 2: 2},
-        8
+        8,
     ],
     # Here the solver has a local minimum to avoid by not choosing the closest point
     "Circles2D3_modified": [
         2,
-        {(1,1): 0, (1,2): 0,
-         (2,1): 4, (2,2): 1,
-         (3,1): 2, (3,2): 4},
+        {(1, 1): 0, (1, 2): 0, (2, 1): 4, (2, 2): 1, (3, 1): 2, (3, 2): 4},
         {1: 1, 2: 1, 3: 1},
         {1: 2, 2: 3, 3: 1},
         {1: 3, 2: 2},
-        8
+        8,
     ],
     # Here's a 3D model.
     "Circles3D4": [
         3,
         {
-            (1, 1): 0, (1, 2): 0, (1, 3): 0,
-            (2, 1): 2, (2, 2): 0, (2, 3): 1,
-            (3, 1): 3, (3, 2): 3, (3, 3): 3,
-            (4, 1): 1, (4, 2): 6, (4, 3): 1,
+            (1, 1): 0,
+            (1, 2): 0,
+            (1, 3): 0,
+            (2, 1): 2,
+            (2, 2): 0,
+            (2, 3): 1,
+            (3, 1): 3,
+            (3, 2): 3,
+            (3, 3): 3,
+            (4, 1): 1,
+            (4, 2): 6,
+            (4, 3): 1,
         },
         {1: 1, 2: 1, 3: 1, 4: 1},
         {1: 1, 2: 3, 3: 2, 4: 1},
         {1: 2, 2: 2, 3: 1},
-        10
-    ]
+        10,
+    ],
 }
 
+
 def build_model(data):
     """Build a model. By default you get Circles2D3, a small instance."""
 
     # Ensure good data was passed
     assert len(data) == 6, "Error processing data"
-    dimension, circle_centers, circle_rvals, circle_penalties, reference_point, upper_bound = data
+    (
+        dimension,
+        circle_centers,
+        circle_rvals,
+        circle_penalties,
+        reference_point,
+        upper_bound,
+    ) = data
 
     assert len(circle_rvals) == len(circle_penalties), "Error processing data"
     assert len(circle_centers) == len(circle_rvals) * dimension, "Error processing data"
@@ -75,37 +87,64 @@ def build_model(data):
     m.circles = RangeSet(len(circle_rvals), doc=f"{len(circle_rvals)} circles")
     m.idx = RangeSet(dimension, doc="n-dimensional indexing set for coordinates")
 
-    m.circ_centers = Param(m.circles, m.idx, initialize=circle_centers, doc="Center points of the circles")
-    m.circ_rvals = Param(m.circles, initialize=circle_rvals, doc="Squared radii of circles")
-    m.ref_point = Param(m.idx, initialize=reference_point, doc="Reference point for distance calculations")
-    m.circ_penalties = Param(m.circles, initialize=circle_penalties, doc="Penalty for being in each circle")
+    m.circ_centers = Param(
+        m.circles, m.idx, initialize=circle_centers, doc="Center points of the circles"
+    )
+    m.circ_rvals = Param(
+        m.circles, initialize=circle_rvals, doc="Squared radii of circles"
+    )
+    m.ref_point = Param(
+        m.idx,
+        initialize=reference_point,
+        doc="Reference point for distance calculations",
+    )
+    m.circ_penalties = Param(
+        m.circles, initialize=circle_penalties, doc="Penalty for being in each circle"
+    )
 
     # Choose a point to minimize objective
-    m.point = Var(m.idx, bounds=(0, upper_bound), doc="Chosen point at which we evaluate objective function")
+    m.point = Var(
+        m.idx,
+        bounds=(0, upper_bound),
+        doc="Chosen point at which we evaluate objective function",
+    )
 
     # Let's set up the "active penalty" this way
-    m.active_penalty = Var(bounds=(0, max(m.circ_penalties[i] for i in m.circles)), doc="Penalty for being in the current circle")
+    m.active_penalty = Var(
+        bounds=(0, max(m.circ_penalties[i] for i in m.circles)),
+        doc="Penalty for being in the current circle",
+    )
 
     # Disjunction: we must be in at least one circle (in fact exactly one since they are disjoint)
     @m.Disjunct(m.circles)
     def circle_disjunct(d, circ):
         m = d.model()
-        d.inside_circle = Constraint(expr=sum((m.point[i] - m.circ_centers[circ, i]) ** 2 for i in m.idx) <= m.circ_rvals[circ])
+        d.inside_circle = Constraint(
+            expr=sum((m.point[i] - m.circ_centers[circ, i]) ** 2 for i in m.idx)
+            <= m.circ_rvals[circ]
+        )
         d.penalty = Constraint(expr=m.active_penalty == m.circ_penalties[circ])
+
     @m.Disjunction()
     def circle_disjunction(m):
         return [m.circle_disjunct[i] for i in m.circles]
 
     # Objective function is Euclidean distance from reference point, plus a penalty depending on which circle we are in
-    m.cost = Objective(expr=sum((m.point[i] - m.ref_point[i]) ** 2 for i in m.idx) + m.active_penalty, doc="Distance from reference point plus penalty")   
+    m.cost = Objective(
+        expr=sum((m.point[i] - m.ref_point[i]) ** 2 for i in m.idx) + m.active_penalty,
+        doc="Distance from reference point plus penalty",
+    )
 
     return m
 
+
 def draw_model(m, title=None):
     """Draw a model using matplotlib to illustrate what's going on. Pass 'title' arg to give chart a title"""
 
     if len(m.idx) != 2:
-        print(f"Unable to draw: this drawing code only supports 2D models, but a {len(m.idx)}-dimensional model was passed.")
+        print(
+            f"Unable to draw: this drawing code only supports 2D models, but a {len(m.idx)}-dimensional model was passed."
+        )
         return
 
     # matplotlib setup
@@ -132,10 +171,7 @@ def draw_model(m, title=None):
 
         ax.add_patch(
             mpl.patches.Circle(
-                (x, y),
-                radius=r,
-                facecolor="#0d98e6",
-                edgecolor="#000000",
+                (x, y), radius=r, facecolor="#0d98e6", edgecolor="#000000"
             )
         )
         # the alignment appears to work by magic
@@ -181,7 +217,6 @@ def draw_model(m, title=None):
     plt.show()
 
 
-
 if __name__ == "__main__":
     from pyomo.environ import SolverFactory
     from pyomo.core.base import TransformationFactory
@@ -189,13 +224,20 @@ def draw_model(m, title=None):
     # Set up a solver, for example scip and bigm works
     solver = SolverFactory("scip")
     transformer = TransformationFactory("gdp.bigm")
-    
 
     for key in circles_model_examples:
         model = build_model(circles_model_examples[key])
         print(f"Solving model {key}")
         transformer.apply_to(model)
         solver.solve(model)
-        pt_string ="(" + "".join(str(value(model.point[i])) + (", " if i != len(model.idx) else "") for i in model.idx) + ")"
+        pt_string = (
+            "("
+            + "".join(
+                str(value(model.point[i])) + (", " if i != len(model.idx) else "")
+                for i in model.idx
+            )
+            + ")"
+        )
         print(f"Optimal value found: {model.cost()} with chosen point {pt_string}")
-        draw_model(model, title=key)
\ No newline at end of file
+        draw_model(model, title=key)
+        print()

From 95c75e37025b28c9ed90ea04b88317cbfb82ab9e Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Tue, 1 Aug 2023 13:13:11 -0600
Subject: [PATCH 09/15] minor readability enhancements on circles example

---
 examples/gdp/circles/circles.py | 68 ++++++++++++++++-----------------
 1 file changed, 34 insertions(+), 34 deletions(-)

diff --git a/examples/gdp/circles/circles.py b/examples/gdp/circles/circles.py
index 74cc74e7957..dfc92357dcb 100644
--- a/examples/gdp/circles/circles.py
+++ b/examples/gdp/circles/circles.py
@@ -18,31 +18,32 @@
 
 # Circles2D3 is the original example described in the papers. Ruiz and Grossman
 # (2012) discuss two more example instances, Circles2D36 and Circles3D36, but I
-# couldn't find them.
+# couldn't find data for them, so I made a couple.
 # format: dimension, circle_centers, circle_rvals, circle_penalties, reference_point, upper_bound
 circles_model_examples = {
-    "Circles2D3": [
-        2,
-        # {1: [0, 0], 2: [4, 1], 3: [2, 4]},
-        {(1, 1): 0, (1, 2): 0, (2, 1): 4, (2, 2): 1, (3, 1): 2, (3, 2): 4},
-        {1: 1, 2: 1, 3: 1},
-        {1: 2, 2: 1, 3: 3},
-        {1: 3, 2: 2},
-        8,
-    ],
+    "Circles2D3": {
+        'dimension': 2,
+        # written as {(circle number, coordinate index): value}, for example below
+        # means {circle 1: (0, 0), circle 2: (4, 1), circle 3: (2, 4)},
+        'circle_centers': {(1, 1): 0, (1, 2): 0, (2, 1): 4, (2, 2): 1, (3, 1): 2, (3, 2): 4},
+        'circle_rvals': {1: 1, 2: 1, 3: 1},
+        'circle_penalties': {1: 2, 2: 1, 3: 3},
+        'reference_point': {1: 3, 2: 2},
+        'upper_bound': 8,
+    },
     # Here the solver has a local minimum to avoid by not choosing the closest point
-    "Circles2D3_modified": [
-        2,
-        {(1, 1): 0, (1, 2): 0, (2, 1): 4, (2, 2): 1, (3, 1): 2, (3, 2): 4},
-        {1: 1, 2: 1, 3: 1},
-        {1: 2, 2: 3, 3: 1},
-        {1: 3, 2: 2},
-        8,
-    ],
+    "Circles2D3_modified": {
+        'dimension': 2,
+        'circle_centers': {(1, 1): 0, (1, 2): 0, (2, 1): 4, (2, 2): 1, (3, 1): 2, (3, 2): 4},
+        'circle_rvals': {1: 1, 2: 1, 3: 1},
+        'circle_penalties': {1: 2, 2: 3, 3: 1},
+        'reference_point': {1: 3, 2: 2},
+        'upper_bound': 8,
+    },
     # Here's a 3D model.
-    "Circles3D4": [
-        3,
-        {
+    "Circles3D4": {
+        'dimension': 3,
+        'circle_centers': {
             (1, 1): 0,
             (1, 2): 0,
             (1, 3): 0,
@@ -56,11 +57,11 @@
             (4, 2): 6,
             (4, 3): 1,
         },
-        {1: 1, 2: 1, 3: 1, 4: 1},
-        {1: 1, 2: 3, 3: 2, 4: 1},
-        {1: 2, 2: 2, 3: 1},
-        10,
-    ],
+        'circle_rvals': {1: 1, 2: 1, 3: 1, 4: 1},
+        'circle_penalties': {1: 1, 2: 3, 3: 2, 4: 1},
+        'reference_point': {1: 2, 2: 2, 3: 1},
+        'upper_bound': 10,
+    },
 }
 
 
@@ -69,14 +70,13 @@ def build_model(data):
 
     # Ensure good data was passed
     assert len(data) == 6, "Error processing data"
-    (
-        dimension,
-        circle_centers,
-        circle_rvals,
-        circle_penalties,
-        reference_point,
-        upper_bound,
-    ) = data
+
+    dimension = data['dimension']
+    circle_centers = data['circle_centers']
+    circle_rvals = data['circle_rvals']
+    circle_penalties = data['circle_penalties']
+    reference_point = data['reference_point']
+    upper_bound = data['upper_bound']
 
     assert len(circle_rvals) == len(circle_penalties), "Error processing data"
     assert len(circle_centers) == len(circle_rvals) * dimension, "Error processing data"

From 61d42023d1e2058e7c819be42f66e43e56fe7e02 Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Tue, 1 Aug 2023 13:14:31 -0600
Subject: [PATCH 10/15] apply black

---
 examples/gdp/circles/circles.py | 18 ++++++++++++++++--
 1 file changed, 16 insertions(+), 2 deletions(-)

diff --git a/examples/gdp/circles/circles.py b/examples/gdp/circles/circles.py
index dfc92357dcb..d0a6a4b1a09 100644
--- a/examples/gdp/circles/circles.py
+++ b/examples/gdp/circles/circles.py
@@ -25,7 +25,14 @@
         'dimension': 2,
         # written as {(circle number, coordinate index): value}, for example below
         # means {circle 1: (0, 0), circle 2: (4, 1), circle 3: (2, 4)},
-        'circle_centers': {(1, 1): 0, (1, 2): 0, (2, 1): 4, (2, 2): 1, (3, 1): 2, (3, 2): 4},
+        'circle_centers': {
+            (1, 1): 0,
+            (1, 2): 0,
+            (2, 1): 4,
+            (2, 2): 1,
+            (3, 1): 2,
+            (3, 2): 4,
+        },
         'circle_rvals': {1: 1, 2: 1, 3: 1},
         'circle_penalties': {1: 2, 2: 1, 3: 3},
         'reference_point': {1: 3, 2: 2},
@@ -34,7 +41,14 @@
     # Here the solver has a local minimum to avoid by not choosing the closest point
     "Circles2D3_modified": {
         'dimension': 2,
-        'circle_centers': {(1, 1): 0, (1, 2): 0, (2, 1): 4, (2, 2): 1, (3, 1): 2, (3, 2): 4},
+        'circle_centers': {
+            (1, 1): 0,
+            (1, 2): 0,
+            (2, 1): 4,
+            (2, 2): 1,
+            (3, 1): 2,
+            (3, 2): 4,
+        },
         'circle_rvals': {1: 1, 2: 1, 3: 1},
         'circle_penalties': {1: 2, 2: 3, 3: 1},
         'reference_point': {1: 3, 2: 2},

From aa799bf4f02d4cc11ea0acabd0826b53f3e3e3f0 Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Tue, 1 Aug 2023 13:29:29 -0600
Subject: [PATCH 11/15] minor improvements for farm layout example

---
 examples/gdp/farm_layout/farm_layout.py | 86 +++++++++++++------------
 1 file changed, 46 insertions(+), 40 deletions(-)

diff --git a/examples/gdp/farm_layout/farm_layout.py b/examples/gdp/farm_layout/farm_layout.py
index 0fc071df90a..3a2ed869793 100644
--- a/examples/gdp/farm_layout/farm_layout.py
+++ b/examples/gdp/farm_layout/farm_layout.py
@@ -22,46 +22,46 @@
 )
 
 # Format: areas, length lower bounds, width lower bounds, length overall upper bound, width overall upper bound
-# Examples are from Sawaya (2006)
+# Examples are from Sawaya (2006), except the labelled alternate ones.
 farm_layout_model_examples = {
-    "FLay02": [{1: 40, 2: 50}, {1: 1, 2: 1}, {1: 1, 2: 1}, 30, 30],
-    "FLay03": [{1: 40, 2: 50, 3: 60}, {1: 1, 2: 1, 3: 1}, {1: 1, 2: 1, 3: 1}, 30, 30],
-    "FLay04": [
-        {1: 40, 2: 50, 3: 60, 4: 35},
-        {1: 3, 2: 3, 3: 3, 4: 3},
-        {1: 3, 2: 3, 3: 3, 4: 3},
-        100,
-        100,
-    ],
-    "FLay05": [
-        {1: 40, 2: 50, 3: 60, 4: 35, 5: 75},
-        {1: 1, 2: 1, 3: 1, 4: 1, 5: 1},
-        {1: 1, 2: 1, 3: 1, 4: 1, 5: 1},
-        30,
-        30,
-    ],
-    "FLay06": [
-        {1: 40, 2: 50, 3: 60, 4: 35, 5: 75, 6: 20},
-        {1: 1, 2: 1, 3: 1, 4: 1, 5: 1, 6: 1},
-        {1: 1, 2: 1, 3: 1, 4: 1, 5: 1, 6: 1},
-        30,
-        30,
-    ],
+    "FLay02": {'areas': {1: 40, 2: 50}, 'length_lbs': {1: 1, 2: 1}, 'width_lbs': {1: 1, 2: 1}, 'length_upper_overall': 30, 'width_upper_overall': 30},
+    "FLay03": {'areas': {1: 40, 2: 50, 3: 60}, 'length_lbs': {1: 1, 2: 1, 3: 1}, 'width_lbs': {1: 1, 2: 1, 3: 1}, 'length_upper_overall': 30, 'width_upper_overall': 30},
+    "FLay04": {
+        'areas': {1: 40, 2: 50, 3: 60, 4: 35},
+        'length_lbs': {1: 3, 2: 3, 3: 3, 4: 3},
+        'width_lbs': {1: 3, 2: 3, 3: 3, 4: 3},
+        'length_upper_overall': 100,
+        'width_upper_overall': 100,
+    },
+    "FLay05": {
+        'areas': {1: 40, 2: 50, 3: 60, 4: 35, 5: 75},
+        'length_lbs': {1: 1, 2: 1, 3: 1, 4: 1, 5: 1},
+        'width_lbs': {1: 1, 2: 1, 3: 1, 4: 1, 5: 1},
+        'length_upper_overall': 30,
+        'width_upper_overall': 30,
+    },
+    "FLay06": {
+        'areas': {1: 40, 2: 50, 3: 60, 4: 35, 5: 75, 6: 20},
+        'length_lbs': {1: 1, 2: 1, 3: 1, 4: 1, 5: 1, 6: 1},
+        'width_lbs': {1: 1, 2: 1, 3: 1, 4: 1, 5: 1, 6: 1},
+        'length_upper_overall': 30,
+        'width_upper_overall': 30,
+    },
     # An example where the dimension constraints are strict enough that it cannot fit a square
-    "FLay03_alt_1": [
-        {1: 40, 2: 50, 3: 60},
-        {1: 5, 2: 5, 3: 5},
-        {1: 5, 2: 5, 3: 7},
-        30,
-        30,
-    ],
-    "FLay03_alt_2": [
-        {1: 40, 2: 50, 3: 60},
-        {1: 1, 2: 1, 3: 1},
-        {1: 1, 2: 1, 3: 1},
-        8,
-        30,
-    ],
+    "FLay03_alt_1": {
+        'areas': {1: 40, 2: 50, 3: 60},
+        'length_lbs': {1: 5, 2: 5, 3: 5},
+        'width_lbs': {1: 5, 2: 5, 3: 7},
+        'length_upper_overall': 30,
+        'width_upper_overall': 30,
+    },
+    "FLay03_alt_2": {
+        'areas': {1: 40, 2: 50, 3: 60},
+        'length_lbs': {1: 1, 2: 1, 3: 1},
+        'width_lbs': {1: 1, 2: 1, 3: 1},
+        'length_upper_overall': 8,
+        'width_upper_overall': 30,
+    },
 }
 
 
@@ -71,8 +71,14 @@ def build_model(params=farm_layout_model_examples["FLay03"]):
     # Ensure good data was passed
     assert (
         len(params) == 5
-    ), "Params should be in the format: areas, length_lbs, width_lbs, length_upper_overall, width_upper, overall"
-    areas, length_lbs, width_lbs, length_upper_overall, width_upper_overall = params
+    ), "Params should be in the format: areas, length_lbs, width_lbs, length_upper_overall, width_upper_overall"
+
+    areas = params['areas']
+    length_lbs = params['length_lbs']
+    width_lbs = params['width_lbs']
+    length_upper_overall = params['length_upper_overall']
+    width_upper_overall = params['width_upper_overall']
+
     for a in areas:
         assert a > 0, "Area must be positive"
 

From 8fe5d8b99e548faf74dfb6c829acebdde764c05d Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Tue, 1 Aug 2023 13:30:39 -0600
Subject: [PATCH 12/15] apply black

---
 examples/gdp/farm_layout/farm_layout.py | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/examples/gdp/farm_layout/farm_layout.py b/examples/gdp/farm_layout/farm_layout.py
index 3a2ed869793..768b127b0da 100644
--- a/examples/gdp/farm_layout/farm_layout.py
+++ b/examples/gdp/farm_layout/farm_layout.py
@@ -24,8 +24,20 @@
 # Format: areas, length lower bounds, width lower bounds, length overall upper bound, width overall upper bound
 # Examples are from Sawaya (2006), except the labelled alternate ones.
 farm_layout_model_examples = {
-    "FLay02": {'areas': {1: 40, 2: 50}, 'length_lbs': {1: 1, 2: 1}, 'width_lbs': {1: 1, 2: 1}, 'length_upper_overall': 30, 'width_upper_overall': 30},
-    "FLay03": {'areas': {1: 40, 2: 50, 3: 60}, 'length_lbs': {1: 1, 2: 1, 3: 1}, 'width_lbs': {1: 1, 2: 1, 3: 1}, 'length_upper_overall': 30, 'width_upper_overall': 30},
+    "FLay02": {
+        'areas': {1: 40, 2: 50},
+        'length_lbs': {1: 1, 2: 1},
+        'width_lbs': {1: 1, 2: 1},
+        'length_upper_overall': 30,
+        'width_upper_overall': 30,
+    },
+    "FLay03": {
+        'areas': {1: 40, 2: 50, 3: 60},
+        'length_lbs': {1: 1, 2: 1, 3: 1},
+        'width_lbs': {1: 1, 2: 1, 3: 1},
+        'length_upper_overall': 30,
+        'width_upper_overall': 30,
+    },
     "FLay04": {
         'areas': {1: 40, 2: 50, 3: 60, 4: 35},
         'length_lbs': {1: 3, 2: 3, 3: 3, 4: 3},

From c4204bcb03a48f3c45d242abd0a39b8544ed0c56 Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Wed, 23 Aug 2023 13:07:09 -0600
Subject: [PATCH 13/15] circles model fixes

---
 examples/gdp/circles/circles.py | 13 +++++--------
 1 file changed, 5 insertions(+), 8 deletions(-)

diff --git a/examples/gdp/circles/circles.py b/examples/gdp/circles/circles.py
index d0a6a4b1a09..ae905998403 100644
--- a/examples/gdp/circles/circles.py
+++ b/examples/gdp/circles/circles.py
@@ -4,8 +4,6 @@
 disjoint hyperspheres.
 """
 
-from __future__ import division
-
 from pyomo.environ import (
     ConcreteModel,
     Objective,
@@ -15,10 +13,13 @@
     Constraint,
     value,
 )
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+
 
 # Circles2D3 is the original example described in the papers. Ruiz and Grossman
 # (2012) discuss two more example instances, Circles2D36 and Circles3D36, but I
-# couldn't find data for them, so I made a couple.
+# (S. Davis 08/2023) couldn't find data for them, so I made a couple.
 # format: dimension, circle_centers, circle_rvals, circle_penalties, reference_point, upper_bound
 circles_model_examples = {
     "Circles2D3": {
@@ -79,7 +80,7 @@
 }
 
 
-def build_model(data):
+def build_model(data=circles_model_examples["Circles2D3"]):
     """Build a model. By default you get Circles2D3, a small instance."""
 
     # Ensure good data was passed
@@ -161,10 +162,6 @@ def draw_model(m, title=None):
         )
         return
 
-    # matplotlib setup
-    import matplotlib as mpl
-    import matplotlib.pyplot as plt
-
     fig, ax = plt.subplots(1, 1, figsize=(10, 10))
     ax.set_xlabel("x")
     ax.set_ylabel("y")

From b13c939811bb36dccf264d3d7f41b47656d6d856 Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Wed, 23 Aug 2023 13:07:36 -0600
Subject: [PATCH 14/15] farm layout model fixes

---
 examples/gdp/farm_layout/farm_layout.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/examples/gdp/farm_layout/farm_layout.py b/examples/gdp/farm_layout/farm_layout.py
index 768b127b0da..f92b3faf376 100644
--- a/examples/gdp/farm_layout/farm_layout.py
+++ b/examples/gdp/farm_layout/farm_layout.py
@@ -8,8 +8,6 @@
 decided that they are on the bottom-left and adapted the disjunction constraints to match.
 """
 
-from __future__ import division
-
 from pyomo.environ import (
     ConcreteModel,
     Objective,
@@ -129,10 +127,10 @@ def build_model(params=farm_layout_model_examples["FLay03"]):
     )
 
     m.plot_x = Var(
-        m.plots, bounds=(0, m.overall_length_ub), doc="x-coordinate of plot center"
+        m.plots, bounds=(0, m.overall_length_ub), doc="x-coordinate of plot bottom-left corner"
     )
     m.plot_y = Var(
-        m.plots, bounds=(0, m.overall_width_ub), doc="y-coordinate of plot center"
+        m.plots, bounds=(0, m.overall_width_ub), doc="y-coordinate of plot bottom-left corner"
     )
     m.plot_length = Var(
         m.plots, bounds=(0, m.overall_length_ub), doc="Length of this plot"

From 45fbf9452e2822e91f58297822d12ab67264b5c9 Mon Sep 17 00:00:00 2001
From: Soren Davis <sorenadavis@gmail.com>
Date: Wed, 23 Aug 2023 13:21:03 -0600
Subject: [PATCH 15/15] apply black

---
 examples/gdp/farm_layout/farm_layout.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/examples/gdp/farm_layout/farm_layout.py b/examples/gdp/farm_layout/farm_layout.py
index f92b3faf376..411e2de3242 100644
--- a/examples/gdp/farm_layout/farm_layout.py
+++ b/examples/gdp/farm_layout/farm_layout.py
@@ -127,10 +127,14 @@ def build_model(params=farm_layout_model_examples["FLay03"]):
     )
 
     m.plot_x = Var(
-        m.plots, bounds=(0, m.overall_length_ub), doc="x-coordinate of plot bottom-left corner"
+        m.plots,
+        bounds=(0, m.overall_length_ub),
+        doc="x-coordinate of plot bottom-left corner",
     )
     m.plot_y = Var(
-        m.plots, bounds=(0, m.overall_width_ub), doc="y-coordinate of plot bottom-left corner"
+        m.plots,
+        bounds=(0, m.overall_width_ub),
+        doc="y-coordinate of plot bottom-left corner",
     )
     m.plot_length = Var(
         m.plots, bounds=(0, m.overall_length_ub), doc="Length of this plot"