diff --git a/examples/gdp/constrained_layout/cons_layout_model.py b/examples/gdp/constrained_layout/cons_layout_model.py index f9110ea5fdf..10595db4c22 100644 --- a/examples/gdp/constrained_layout/cons_layout_model.py +++ b/examples/gdp/constrained_layout/cons_layout_model.py @@ -11,29 +11,116 @@ """ from __future__ import division -from pyomo.environ import ConcreteModel, Objective, Param, RangeSet, Set, Var +from pyomo.environ import ConcreteModel, Objective, Param, RangeSet, Set, Var, value +# Constrained layout model examples. These are from Nicolas Sawaya (2006). +# Format: rect_lengths, rect_heights, circ_xvals, circ_yvals, circ_rvals (as dicts), +# sep_penalty_matrix (as nested array) +# Note that only the strict upper triangle of sep_penalty_matrix is used +constrained_layout_model_examples = { + "CLay0203": { + 'rect_lengths': {1: 5, 2: 7, 3: 3}, + 'rect_heights': {1: 6, 2: 5, 3: 3}, + 'circ_xvals': {1: 15, 2: 50}, + 'circ_yvals': {1: 10, 2: 80}, + 'circ_rvals': {1: 6, 2: 5}, + 'sep_penalty_matrix': [[0, 300, 240], [0, 0, 100]], + }, + "CLay0204": { + 'rect_lengths': {1: 5, 2: 7, 3: 3, 4: 2}, + 'rect_heights': {1: 6, 2: 5, 3: 3, 4: 3}, + 'circ_xvals': {1: 15, 2: 50}, + 'circ_yvals': {1: 10, 2: 80}, + 'circ_rvals': {1: 6, 2: 10}, + 'sep_penalty_matrix': [[0, 300, 240, 210], [0, 0, 100, 150], [0, 0, 0, 120]], + }, + "CLay0205": { + 'rect_lengths': {1: 5, 2: 7, 3: 3, 4: 2, 5: 9}, + 'rect_heights': {1: 6, 2: 5, 3: 3, 4: 3, 5: 7}, + 'circ_xvals': {1: 15, 2: 50}, + 'circ_yvals': {1: 10, 2: 80}, + 'circ_rvals': {1: 6, 2: 10}, + 'sep_penalty_matrix': [ + [0, 300, 240, 210, 50], + [0, 0, 100, 150, 30], + [0, 0, 0, 120, 25], + [0, 0, 0, 0, 60], + ], + }, + "CLay0303": { + 'rect_lengths': {1: 5, 2: 7, 3: 3}, + 'rect_heights': {1: 6, 2: 5, 3: 3}, + 'circ_xvals': {1: 15, 2: 50, 3: 30}, + 'circ_yvals': {1: 10, 2: 80, 3: 50}, + 'circ_rvals': {1: 6, 2: 5, 3: 4}, + 'sep_penalty_matrix': [[0, 300, 240], [0, 0, 100]], + }, + "CLay0304": { + 'rect_lengths': {1: 5, 2: 7, 3: 3, 4: 2}, + 'rect_heights': {1: 6, 2: 5, 3: 3, 4: 3}, + 'circ_xvals': {1: 15, 2: 50, 3: 30}, + 'circ_yvals': {1: 10, 2: 80, 3: 50}, + 'circ_rvals': {1: 6, 2: 5, 3: 4}, + 'sep_penalty_matrix': [[0, 300, 240, 210], [0, 0, 100, 150], [0, 0, 0, 120]], + }, + "CLay0305": { + 'rect_lengths': {1: 5, 2: 7, 3: 3, 4: 2, 5: 9}, + 'rect_heights': {1: 6, 2: 5, 3: 3, 4: 3, 5: 7}, + 'circ_xvals': {1: 15, 2: 50, 3: 30}, + 'circ_yvals': {1: 10, 2: 80, 3: 50}, + 'circ_rvals': {1: 6, 2: 10, 3: 4}, + 'sep_penalty_matrix': [ + [0, 300, 240, 210, 50], + [0, 0, 100, 150, 30], + [0, 0, 0, 120, 25], + [0, 0, 0, 0, 60], + ], + }, +} -def build_constrained_layout_model(): + +def build_constrained_layout_model( + params=constrained_layout_model_examples['CLay0203'], metric="l1" +): """Build the model.""" + + # Ensure the caller passed good data + assert len(params) == 6 + + # Get all the parameters out + rect_lengths = params['rect_lengths'] + rect_heights = params['rect_heights'] + circ_xvals = params['circ_xvals'] + circ_yvals = params['circ_yvals'] + circ_rvals = params['circ_rvals'] + sep_penalty_matrix = params['sep_penalty_matrix'] + + assert len(rect_lengths) == len( + rect_heights + ), "There should be the same number of rectangle lengths and heights." + assert ( + len(circ_xvals) == len(circ_yvals) == len(circ_rvals) + ), "There should be the same number of circle x values, y values, and radii" + for row in sep_penalty_matrix: + assert len(row) == len( + sep_penalty_matrix[0] + ), "Matrix rows should have the same length" + assert metric in ["l1", "l2"], 'Metric options are "l1" and "l2"' + m = ConcreteModel(name="2-D constrained layout") - m.rectangles = RangeSet(3, doc="Three rectangles") - m.circles = RangeSet(2, doc="Two circles") + m.rectangles = RangeSet(len(rect_lengths), doc=f"{len(rect_lengths)} rectangles") + m.circles = RangeSet(len(circ_xvals), doc=f"{len(circ_xvals)} circles") - m.rect_length = Param( - m.rectangles, initialize={1: 5, 2: 7, 3: 3}, doc="Rectangle length" - ) - m.rect_height = Param( - m.rectangles, initialize={1: 6, 2: 5, 3: 3}, doc="Rectangle height" - ) + m.rect_length = Param(m.rectangles, initialize=rect_lengths, doc="Rectangle length") + m.rect_height = Param(m.rectangles, initialize=rect_heights, doc="Rectangle height") m.circle_x = Param( - m.circles, initialize={1: 15, 2: 50}, doc="x-coordinate of circle center" + m.circles, initialize=circ_xvals, doc="x-coordinate of circle center" ) m.circle_y = Param( - m.circles, initialize={1: 10, 2: 80}, doc="y-coordinate of circle center" + m.circles, initialize=circ_yvals, doc="y-coordinate of circle center" ) - m.circle_r = Param(m.circles, initialize={1: 6, 2: 5}, doc="radius of circle") + m.circle_r = Param(m.circles, initialize=circ_rvals, doc="radius of circle") @m.Param(m.rectangles, doc="Minimum feasible x value for rectangle") def x_min(m, rect): @@ -63,16 +150,16 @@ def y_max(m, rect): for circ in m.circles ) - m.ordered_rect_pairs = Set( - initialize=m.rectangles * m.rectangles, filter=lambda _, r1, r2: r1 != r2 - ) m.rect_pairs = Set( - initialize=[(r1, r2) for r1, r2 in m.ordered_rect_pairs if r1 < r2] + initialize=m.rectangles * m.rectangles, filter=lambda _, r1, r2: r1 < r2 ) m.rect_sep_penalty = Param( m.rect_pairs, - initialize={(1, 2): 300, (1, 3): 240, (2, 3): 100}, + # 0-based vs 1-based indices... + initialize={ + (r1, r2): sep_penalty_matrix[r1 - 1][r2 - 1] for r1, r2 in m.rect_pairs + }, doc="Penalty for separation distance between two rectangles.", ) @@ -91,20 +178,40 @@ def y_bounds(m, rect): m.dist_x = Var(m.rect_pairs, doc="x-axis separation between rectangle pair") m.dist_y = Var(m.rect_pairs, doc="y-axis separation between rectangle pair") - m.min_dist_cost = Objective( - expr=sum( - m.rect_sep_penalty[r1, r2] * (m.dist_x[r1, r2] + m.dist_y[r1, r2]) - for (r1, r2) in m.rect_pairs + if metric == "l2": + m.min_dist_cost = Objective( + expr=sum( + m.rect_sep_penalty[r1, r2] + * (m.dist_x[r1, r2] ** 2 + m.dist_y[r1, r2] ** 2) ** 0.5 + for (r1, r2) in m.rect_pairs + ) ) - ) + # l1 distance used in the paper + else: + m.min_dist_cost = Objective( + expr=sum( + m.rect_sep_penalty[r1, r2] * (m.dist_x[r1, r2] + m.dist_y[r1, r2]) + for (r1, r2) in m.rect_pairs + ) + ) + + # Ensure the dist_x and dist_y are greater than the positive and negative + # signed distances. + @m.Constraint(m.rect_pairs, doc="x-distance between rectangles") + def dist_x_defn_1(m, r1, r2): + return m.dist_x[r1, r2] >= m.rect_x[r2] - m.rect_x[r1] + + @m.Constraint(m.rect_pairs, doc="x-distance between rectangles") + def dist_x_defn_2(m, r1, r2): + return m.dist_x[r1, r2] >= m.rect_x[r1] - m.rect_x[r2] - @m.Constraint(m.ordered_rect_pairs, doc="x-distance between rectangles") - def dist_x_defn(m, r1, r2): - return m.dist_x[tuple(sorted([r1, r2]))] >= m.rect_x[r2] - m.rect_x[r1] + @m.Constraint(m.rect_pairs, doc="y-distance between rectangles") + def dist_y_defn_1(m, r1, r2): + return m.dist_y[r1, r2] >= m.rect_y[r2] - m.rect_y[r1] - @m.Constraint(m.ordered_rect_pairs, doc="y-distance between rectangles") - def dist_y_defn(m, r1, r2): - return m.dist_y[tuple(sorted([r1, r2]))] >= m.rect_y[r2] - m.rect_y[r1] + @m.Constraint(m.rect_pairs, doc="y-distance between rectangles") + def dist_y_defn_2(m, r1, r2): + return m.dist_y[r1, r2] >= m.rect_y[r1] - m.rect_y[r2] @m.Disjunction( m.rect_pairs, @@ -148,3 +255,91 @@ def rectangle_in_circle(m, r): ] return m + + +def draw_model(m, title=None): + """Draw a model using matplotlib to illustrate what's going on. Pass 'title' argument 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") + + # Hardcode these bounds since I'm not sure the best way to do it automatically + ax.set_xlim([0, 100]) + ax.set_ylim([0, 100]) + + if title is not None: + plt.title(title) + + for c in m.circles: + print( + f"drawing circle {c}: x={m.circle_x[c]}, y={m.circle_y[c]}, r={m.circle_r[c]}" + ) + # no idea about colors + ax.add_patch( + mpl.patches.Circle( + (m.circle_x[c], m.circle_y[c]), + radius=m.circle_r[c], + facecolor="#0d98e6", + edgecolor="#000000", + ) + ) + ax.text( + m.circle_x[c], + m.circle_y[c] + m.circle_r[c] + 1.5, + f"Circle {c}", + horizontalalignment="center", + ) + + for r in m.rectangles: + print( + f"drawing rectangle {r}: x={value(m.rect_x[r])}, y={value(m.rect_y[r])} (center), L={m.rect_length[r]}, H={m.rect_height[r]}" + ) + ax.add_patch( + mpl.patches.Rectangle( + ( + value(m.rect_x[r]) - m.rect_length[r] / 2, + value(m.rect_y[r]) - m.rect_height[r] / 2, + ), + m.rect_length[r], + m.rect_height[r], + facecolor="#fbec68", + edgecolor="#000000", + ) + ) + ax.text( + value(m.rect_x[r]), + value(m.rect_y[r]), + f"R{r}", + horizontalalignment="center", + verticalalignment="center", + fontsize=8, + ) + + plt.show() + + +# Solve some example problems and draw some pictures +if __name__ == "__main__": + from pyomo.environ import SolverFactory, TransformationFactory + + # Set up a solver, for example scip and bigm works + solver = SolverFactory("scip") + transformer = TransformationFactory("gdp.bigm") + + # Do all of the possible problems + for d in ["l1", "l2"]: + for key in constrained_layout_model_examples.keys(): + print(f"Solving example problem: {key}") + model = build_constrained_layout_model( + constrained_layout_model_examples[key], metric=d + ) + transformer.apply_to(model) + solver.solve(model) + print(f"Found objective function value: {model.min_dist_cost()}") + draw_model(model, title=(f"{key} ({d} distance)")) + print()