diff --git a/examples/gdp/circles/circles.py b/examples/gdp/circles/circles.py new file mode 100644 index 00000000000..ae905998403 --- /dev/null +++ b/examples/gdp/circles/circles.py @@ -0,0 +1,254 @@ +""" +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 pyomo.environ import ( + ConcreteModel, + Objective, + Param, + RangeSet, + Var, + 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 +# (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": { + '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": { + '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": { + 'dimension': 3, + 'circle_centers': { + (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, + }, + '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, + }, +} + + +def build_model(data=circles_model_examples["Circles2D3"]): + """Build a model. By default you get Circles2D3, a small instance.""" + + # Ensure good data was passed + assert len(data) == 6, "Error processing 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" + 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 + + 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) + print() diff --git a/examples/gdp/farm_layout/farm_layout.py b/examples/gdp/farm_layout/farm_layout.py new file mode 100644 index 00000000000..411e2de3242 --- /dev/null +++ b/examples/gdp/farm_layout/farm_layout.py @@ -0,0 +1,276 @@ +""" +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 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), 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, + }, + "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": { + '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, + }, +} + + +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 = 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" + + 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 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 bottom-left corner", + ) + m.plot_y = Var( + 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" + ) + 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 + @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] <= 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 + # 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]) + + # 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] + + 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""" + + # 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, + ) + ) + + # 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_y[p]()), + m.plot_length[p](), + m.plot_width[p](), + facecolor="#ebdc78", + edgecolor="#000000", + ) + ) + ax.text( + 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", + ) + + plt.show() + + +if __name__ == "__main__": + from pyomo.environ import SolverFactory + from pyomo.core.base import TransformationFactory + + # 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.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()