| High-level model |
|---|
| \[\begin{align} \min\> &{\color{darkred}{\mathit{area}}} = \sum_k \prod_c \left({\color{darkred}{\mathit{rect}}}_{k,c,up}-{\color{darkred}{\mathit{rect}}}_{k,c,lo}\right) \\ & {\color{darkred}\delta}_{i,k} = 1 \implies {\color{darkblue}p}_{i,c} \in [{\color{darkred}{\mathit{rect}}}_{k,c,lo},{\color{darkred}{\mathit{rect}}}_{k,c,up}] & \forall i,k,c\\ & \sum_k{\color{darkred}\delta}_{i,k} = {\color{darkblue}1} & \forall i \\ & {\color{darkred}\delta}_{i,k} \in \{0,1\} \\ & c \in \{x,y\} \end{align}\] |
| Implemented model |
|---|
| \[\begin{align} \min\> &{\color{darkred}{\mathit{area}}} = \sum_k {\color{darkred}d}_{k,x} \cdot{\color{darkred}d}_{k,y} \\ & {\color{darkred}d}_{k,c} = {\color{darkred}{\mathit{rect}}}_{k,c,up}- {\color{darkred}{\mathit{rect}}}_{k,c,lo} & \forall k,c \\ & {\color{darkred}{\mathit{rect}}}_{k,c,lo} \le {\color{darkred}\delta}_{i,k} {\color{darkblue}p}_{i,c} + (1-{\color{darkred}\delta}_{i,k}) {\color{darkblue}{\mathit{box}}}_{c,up} & \forall i,c \\ & {\color{darkred}{\mathit{rect}}}_{k,c,up} \ge {\color{darkred}\delta}_{i,k} {\color{darkblue}p}_{i,c} + (1-{\color{darkred}\delta}_{i,k}){\color{darkblue}{\mathit{box}}}_{c,lo}& \forall i,c\\ & \sum_k{\color{darkred}\delta}_{i,k} = 1 & \forall i \\& {\color{darkred}{\mathit{rect}}}_{k,c,lo},{\color{darkred}{\mathit{rect}}}_{k,c,up} \in [{\color{darkblue}{\mathit{box}}}_{c,lo},{\color{darkblue}{\mathit{box}}}_{c,up}] \\ & {\color{darkred}d}_{k,c} \in [0,{\color{darkblue}{\mathit{box}}}_{c,up}-{\color{darkblue}{\mathit{box}}}_{c,lo}] \\ & {\color{darkred}\delta}_{i,k} \in \{0,1\}\\ &{\color{darkblue}{\mathit{box}}}_{c,lo} = \min_i{\color{darkblue}p}_{i,c} \\ &{\color{darkblue}{\mathit{box}}}_{c,up} = \max_i{\color{darkblue}p}_{i,c}\\ & c \in \{x,y\} \end{align}\] |
The perimeter objective would be: \[\min\>\sum_k 2\left({\color{darkred}d}_{k,x} +{\color{darkred}d}_{k,y} \right)\]We can add to the model a simple symmetry breaker: \[ {\color{darkred}{\mathit{rect}}}_{k,x,lo} \le {\color{darkred}{\mathit{rect}}}_{k+1,x,lo}\] This constraint orders the rectangles by their \(x\)-coordinate. It makes the solution more readable, but more importantly, it reduces symmetry in the model and can increase solver performance, sometimes by a large amount.
- The nonconvex quadratic model is not so easy to solve. You need a global solver for this. The extra symmetry-breaking constraint helps a lot.
- The data is randomly generated. Since we can easily scale this problem, it is a good model for benchmarking and stress testing. I found issues with some solvers.
- The linear MIP model is much easier.
- We could enumerate the possible rectangles. The boundaries always fall on a data point, so this is discrete and can be enumerated. The resulting optimization problem is much larger but much simpler.
Data
Dataset
---- 37 SET k rectangles rectangle1, rectangle2, rectangle3, rectangle4, rectangle5 ---- 37 PARAMETER p points x y point1 4.97311132 9.63862756 point2 8.73922492 8.44273500 point3 1.91226253 4.94526579 point4 2.34553216 3.80983451 point5 6.00887463 1.51476321 point6 5.07716895 4.63323882 point7 3.59064347 1.52625715 point8 3.66418827 8.26045475 point9 0.23893508 0.59261597 point10 2.84759145 1.14012489 point11 8.90454986 3.74340300 point12 2.34069434 1.47696552 point13 0.19119162 0.78401936 point14 9.27971966 5.55039788 point15 6.16361826 3.60477909 point16 5.21416943 4.32016162 point17 7.49929678 8.17376279 point18 0.16037081 3.72949918 point19 6.97091674 7.66895295 point20 8.87457833 8.11673997 point21 7.89650246 8.03634493 point22 0.16576575 3.93244606 point23 9.60138480 7.21520825 point24 8.14734411 8.37501191 point25 7.74541956 5.28258128 point26 6.55663865 3.92384422 point27 4.85091115 4.02642582 point28 4.34623796 1.93472329 point29 1.28959803 4.87601902 point30 7.91163148 6.85615242 point31 3.02621818 0.00063517 point32 4.46278421 8.14820541 point33 2.08022534 3.35953314 point34 2.81600217 1.56499252 point35 6.20283636 8.59944445 point36 6.15258941 3.85719583 point37 7.63343999 7.31506452 point38 0.79832230 0.82942221 point39 0.19238302 2.45380498 point40 3.34088869 3.03838803 point41 8.34902953 3.86320590 point42 0.53110352 5.05611011 point43 1.91904761 2.34489699 point44 9.34705030 7.86066922 point45 9.43453787 1.71763581 point46 1.81723665 2.02565095 point47 5.32342079 5.06474382 point48 2.10786534 6.38173925 point49 2.95917598 5.53252693 point50 0.34180259 8.07512765
Solutions
Run 1: minimize area
---- 113 PARAMETER assign assignment=round(delta)
rectangle1 rectangle2 rectangle3 rectangle4 rectangle5
point1 1
point2 1
point3 1
point4 1
point5 1
point6 1
point7 1
point8 1
point9 1
point10 1
point11 1
point12 1
point13 1
point14 1
point15 1
point16 1
point17 1
point18 1
point19 1
point20 1
point21 1
point22 1
point23 1
point24 1
point25 1
point26 1
point27 1
point28 1
point29 1
point30 1
point31 1
point32 1
point33 1
point34 1
point35 1
point36 1
point37 1
point38 1
point39 1
point40 1
point41 1
point42 1
point43 1
point44 1
point45 1
point46 1
point47 1
point48 1
point49 1
point50 1
count 20 3 4 5 18
We have two large clusters (18 and 20 points) and three small ones (3, 4, and 5 points). I rounded the solution because MIP solvers often return very small values, such as 1e-10. Here we filtered those out. We also added the counts.
If we look at the rectangles, we see:
---- 141 PARAMETER results by cluster INDEX 1 = min area obj min max len points area perimeter rectangle1.x 0.160 3.341 3.181 rectangle1.y 0.001 6.382 6.381 rectangle1. 20.000 20.295 19.123 rectangle2.x 0.342 4.463 4.121 rectangle2.y 8.075 8.260 0.185 rectangle2. 3.000 0.764 8.613 rectangle3.x 3.591 9.435 5.844 rectangle3.y 1.515 1.935 0.420 rectangle3. 4.000 2.454 12.528 rectangle4.x 4.851 5.323 0.473 rectangle4.y 4.026 9.639 5.612 rectangle4. 5.000 2.652 12.169 rectangle5.x 6.153 9.601 3.449 rectangle5.y 3.605 8.599 4.995 rectangle5. 18.000 17.226 16.887 total . 50.000 43.391 69.320
Obviously, the clusters with many points are also the largest (by area or perimeter).
The highlighted number is the one the model tries to minimize.
As a side note: I spent a little bit of time on this table. In general, for more complex models, with more complex results, I think it is always a good idea to think a bit about reporting. Reports, in essence, are the model's user interface (GUI). So they deserve some attention.
A picture of the results is:
This looks quite reasonable.
Run 2: minimize perimeter
The "minimize perimeter" is an easier, linear MIP model. Let's see what happens when we run this version. I expect we will observe a solution close to the one in the previous section.
The optimal assignment is as follows:
---- 133 PARAMETER assign assignment=round(delta)
rectangle1 rectangle2 rectangle3 rectangle4 rectangle5
point1 1
point2 1
point3 1
point4 1
point5 1
point6 1
point7 1
point8 1
point9 1
point10 1
point11 1
point12 1
point13 1
point14 1
point15 1
point16 1
point17 1
point18 1
point19 1
point20 1
point21 1
point22 1
point23 1
point24 1
point25 1
point26 1
point27 1
point28 1
point29 1
point30 1
point31 1
point32 1
point33 1
point34 1
point35 1
point36 1
point37 1
point38 1
point39 1
point40 1
point41 1
point42 1
point43 1
point44 1
point45 1
point46 1
point47 1
point48 1
point49 1
point50 1
count 46 1 1 1 1
The details about the clusters are:
---- 148 PARAMETER results by cluster INDEX 1 = min perim obj min max len points area perimeter rectangle1.x 0.160 9.435 9.274 rectangle1.y 0.784 8.599 7.815 rectangle1. 46.000 72.482 34.179 rectangle2.x 0.239 0.239 rectangle2.y 0.593 0.593 rectangle2. 1.000 rectangle3.x 3.026 3.026 rectangle3.y 0.001 0.001 rectangle3. 1.000 rectangle4.x 4.973 4.973 rectangle4.y 9.639 9.639 rectangle4. 1.000 rectangle5.x 9.601 9.601 rectangle5.y 7.215 7.215 rectangle5. 1.000 total . 50.000 72.482 34.179
This looks very different: we have one big cluster, and then some singletons.
Computation: Symmetry breaker
| nonconvex MIQP model | |
|---|---|
| Without ordering constraint | ... Status: Normal (1) / Global Optimal (1)
... Objective Value: 4.339054584504e+01
... Objective Bound: 4.339054584504e+01
... Optimality Gap: 0.000000000000e+00 (0.00%)
... Solved in 3254.985s, 0 iterations, 11328815 nodes |
| With ordering constraint | ... Status: Normal (1) / Global Optimal (1)
... Objective Value: 4.339054584504e+01
... Objective Bound: 4.339054584504e+01
... Optimality Gap: 0.000000000000e+00 (0.00%)
... Solved in 398.321s, 0 iterations, 1864195 nodes
|
To prove global optimality, adding the symmetry-breaking ordering constraint reduces the solution time by roughly a factor of 8.
The zero-iterations message is obviously very hard to believe. Let's call that a minor bug.
It does not not always work this way. Here are some results for the easier, linear MIP (minimize perimeter):
| linear MIP model | |
|---|---|
| Without ordering constraint | ... Status: Normal (1) / Global Optimal (1)
... Objective Value: 3.417918430000e+01
... Objective Bound: 3.417918430000e+01
... Optimality Gap: 0.000000000000e+00 (0.00%)
... Solved in 0.818s, 108681 iterations, 417 nodes |
| With ordering constraint | ... Status: Normal (1) / Global Optimal (1) ... Objective Value: 3.417918430000e+01 ... Objective Bound: 3.417918430000e+01 ... Optimality Gap: 0.000000000000e+00 (0.00%) ... Solved in 4.926s, 829487 iterations, 8251 nodes |
Run 3: enumeration + set covering approach
In this section, we provide an alternative approach for both the "min area" and "min perimeter" models.
We'll generate all possible rectangles given our set of points \({\color{darkblue}p}_{i,c}\). We can do this efficiently because of the special property that points coincide with the four sides of the rectangles: we enumerate combinations of points that constitute a rectangle. We'll use Python for this and ask Gemini to help generate the code. This actually worked surprisingly well. The generated code isn't trivial.
For our dataset with \(n=50\) points, we generate \(51,587\) rectangles! This is not small.
We can simply pre-calculate the area (or perimeter) for all rectangles, and now we can solve the very large but very easy linear MIP model:
| Linear MIP Model |
|---|
| \[\begin{align}\min& \sum_r {\color{darkblue}{\mathit{area}}}_r \cdot {\color{darkred}{\mathit{select}}}_r \\ & \sum_{r|{\color{darkblue}{\mathit{contain}}}_{i,r}} {\color{darkred}{\mathit{select}}}_r \ge 1 & \forall i \\ & \sum_r {\color{darkred}{\mathit{select}}}_r = \color{darkblue}k \\ & {\color{darkred}{\mathit{select}}}_r \in \{0,1\} \end{align} \] |
Here the set \( {\color{darkblue}{\mathit{contain}}}_{i,r} \) indicates if point \(i\) is contained in rectangle \(r\). You can also implement this with a boolean parameter instead of a set: \[\sum_r {\color{darkblue}{\mathit{contain}}}_{i,r} \cdot {\color{darkred}{\mathit{select}}}_r \ge 1\]
Of course, we can replace the area-based objective with the sum of the perimeters.
If we run this, we see:
| Linear set covering models | |
|---|---|
| Min total area | Best solution : 43.390545845
Best bound : 43.390545845
Best gap : 0.0000%
Solve time : 1.88
Solve node : 1
MIP status : solved
Solution status : integer optimal (relative gap limit 0.0001)
|
| Min total perimeter | Best solution : 34.179184300
Best bound : 34.179184300
Best gap : 0.0000%
Solve time : 8.90
Solve node : 1
MIP status : solved
Solution status : integer optimal (relative gap limit 0.0001)
|
The objectives match what we found before. Only this is way, way faster for the expensive solve. This is a good example that shows the number of integer variables can be a poor predictor of performance.
To summarize, the solution times for the min area problem are:
| Min area solution times (s) | |
|---|---|
| Nonconvex MIQCP | 3254.99 |
| Nonconvex MIQCP + symmetry breaker | 398.32 |
| Enumerate rectangles | 2.59 |
| Set covering MIP | 1.88 |
For the set covering approach, we should add the time to generate the rectangles.
Conclusion
As an aside, we also showed the power of a set-covering approach for these problems. Compared with the rather difficult non-convex MIQCP model, this yields large performance gains. Looking at a problem differently can help.
I found this exercise interesting. The models are nice, clean, and readable. And they make a good point: what people often think of as models that deliver very similar solutions sometimes don't. Sometimes, very different modeling techniques can lead to much better-performing models.
Also, because we use a random dataset and can scale it to any size, this could create good candidates for global benchmark and quality-assurance problems. We found some issues with solvers.
References
- Minimum enclosing rectangle for \(k\) out of \(n\) points, https://yetanothermathprogrammingconsultant.blogspot.com/2026/09/minimum-enclosing-rectangle-for-k-out.html
- K-medians clustering, https://en.wikipedia.org/wiki/K-medians_clustering
Appendix
GAMS Models Run 1,2
$onText
cluster points using axis-aligned rectangles Two implementations: minimize area (nonconvex MIQCP) minimize perimeter (linear MIP)
$offText
*--------------------------------------------------------------- * Settings *---------------------------------------------------------------
option seed=12345;
option mip = xpress; option miqcp = xpress;
*--------------------------------------------------------------- * data *---------------------------------------------------------------
scalars maxsize 'size of the outer box' / 10 / ;
sets i 'points' /point1*point50/ k 'rectangles' /rectangle1*rectangle5/ c 'coordinates' /x,y/ lu 'lower/upper values' /min,max/ ;
parameter p(i,c) 'points'; p(i,c) = uniform(0,maxsize);
option k:0,p:8; display k,p;
*--------------------------------------------------------------- * derived data *---------------------------------------------------------------
* bounding box for data points parameter box(c,lu) 'min/max values'; box(c,'min') = smin(i, p(i,c)); box(c,'max') = smax(i, p(i,c));
display box;
*--------------------------------------------------------------- * non-convex quadratic and linear model *---------------------------------------------------------------
binary variables delta(i,k) 'assign point to rectangle';
positive variables rect(k,c,lu) 'rectangle coordinates' d(k,c) 'width and length of rectangle' ;
variables totalArea 'objective (nonconvex quadratic)' totalPerimeter 'objective (linear)' ;
* bounds rect.lo(k,c,lu) = box(c,'min'); rect.up(k,c,lu) = box(c,'max'); d.up(k,c) = box(c,'max')-box(c,'min');
Equations obj1 'nonconvex quadratic: min sum area' obj2 'linear: min sum perimeter' inside_min(i,k,c) 'point i is inside rectangle k' inside_max(i,k,c) 'point i is inside rectangle k' inside_one(i) 'each point must be inside one rectangle' calc_d(k,c) 'calc height/width of rectangle' extra(k) 'symmetry breaker' ;
obj1.. totalArea =e= sum(k, prod(c,d(k,c))); obj2.. totalPerimeter =e= sum(k, 2*sum(c,d(k,c)));
inside_min(i,k,c).. rect(k,c,'min') =l= delta(i,k)*p(i,c) + (1-delta(i,k))*box(c,'max'); inside_max(i,k,c).. rect(k,c,'max') =g= delta(i,k)*p(i,c) + (1-delta(i,k))*box(c,'min');
inside_one(i).. sum(k, delta(i,k)) =e= 1;
calc_d(k,c).. d(k,c) =e= rect(k,c,'max') - rect(k,c,'min');
extra(k+1).. rect(k,'x','min') =l= rect(k+1,'x','min');
model mArea 'quadratic' / all - obj2 /; model mPerimeter 'linear' / all - obj1 /;
*--------------------------------------------------------------- * Solve and save results *---------------------------------------------------------------
set dummy /min,max,len,points,count,area,perimeter/;
parameter results(*,*,*,*) 'by cluster' assign(*,k) 'assignment=round(delta)' rect1 'saved version for postprocessing' rect2 'saved version for postprocessing' ;
option assign:0;
solve mArea minimizing totalArea using miqcp; abort$(mArea.solvestat<>%solveStat.normalCompletion%) "Abnormal completion"; assign(i,k) = round(delta.l(i,k)); assign('count',k) = sum(i,assign(i,k)); display delta.l, assign, rect.l, d.l, totalArea.l;
rect1(k,c,lu) = rect.l(k,c,lu);
results('min area obj',k,c,'min') = rect.l(k,c,'min'); results('min area obj',k,c,'max') = rect.l(k,c,'max'); results('min area obj',k,c,'len') = d.l(k,c); results('min area obj',k,'','points') = sum(i$(delta.l(i,k)>0.5),1); results('min area obj',k,'','area') = prod(c,d.l(k,c)); results('min area obj',k,'','perimeter') = 2*sum(c,d.l(k,c)); results('min area obj','total','','points') = sum((i,k)$(delta.l(i,k)>0.5),1); results('min area obj','total','','area') = sum(k,prod(c,d.l(k,c))); results('min area obj','total','','perimeter') = sum(k,2*sum(c,d.l(k,c)));
solve mPerimeter minimizing totalPerimeter using mip; abort$(mPerimeter.solvestat<>%solveStat.normalCompletion%) "Abnormal completion"; assign(i,k) = round(delta.l(i,k)); assign('count',k) = sum(i,assign(i,k)); display delta.l, assign, rect.l, d.l, totalPerimeter.l;
rect2(k,c,lu) = rect.l(k,c,lu);
results('min perim obj',k,c,'min') = rect.l(k,c,'min'); results('min perim obj',k,c,'max') = rect.l(k,c,'max'); results('min perim obj',k,c,'len') = d.l(k,c); results('min perim obj',k,'','points') = sum(i$(delta.l(i,k)>0.5),1); results('min perim obj',k,'','area') = prod(c,d.l(k,c)); results('min perim obj',k,'','perimeter') = 2*sum(c,d.l(k,c)); results('min perim obj','total','','points') = sum((i,k)$(delta.l(i,k)>0.5),1); results('min perim obj','total','','area') = sum(k,prod(c,d.l(k,c))); results('min perim obj','total','','perimeter') = sum(k,2*sum(c,d.l(k,c)));
option results:3; display results; |
GAMS/Python code Run 3
Python code: enumeration of rectangles
# renumerate.py version = "1.0.0" import bisect from collections import defaultdict # this code was generated by gemini def find_all_boundary_rectangles(points_list): """ Finds all valid rectangles having a point on each of its 4 boundaries. Parameters: points_list (list): A list of tuples (id: string, x: float, y: float) Returns: dict: A mapping of rectangle IDs to their bounds and covered point IDs. """ # 1. Store points in a dictionary by their ID points_dict = {pid: {"x": px, "y": py} for pid, px, py in points_list} # 2. Group points by X coordinate for efficient slicing x_to_pts = defaultdict(list) for pid, px, py in points_list: x_to_pts[px].append((py, pid)) # Store as (y, id) to sort by Y easily unique_x = sorted(list(x_to_pts.keys())) rectangles_output = {} rect_counter = 1 # Iterate over all pairs of X boundaries for i in range(len(unique_x)): x_min = unique_x[i] # Pull points incrementally between x_min and x_max current_pts = [] for j in range(i, len(unique_x)): x_max = unique_x[j] current_pts.extend(x_to_pts[x_max]) # Left and right boundary Y arrays for fast constraint checks left_ys = sorted([item[0] for item in x_to_pts[x_min]]) right_ys = sorted([item[0] for item in x_to_pts[x_max]]) # Sort current sweep-line slice by Y-coordinate current_pts.sort(key=lambda item: item[0]) all_unique_ys = sorted(list(set(item[0] for item in current_pts))) start_pt_idx = 0 for y_min in all_unique_ys: # Find smallest left Y >= y_min idx_l = bisect.bisect_left(left_ys, y_min) if idx_l == len(left_ys): continue # Find smallest right Y >= y_min idx_r = bisect.bisect_left(right_ys, y_min) if idx_r == len(right_ys): continue # Minimum height constraint to close left/right boundaries y_start = max(left_ys[idx_l], right_ys[idx_r]) y_max_choices = [y for y in all_unique_ys if y >= y_start] if not y_max_choices: continue while start_pt_idx < len(current_pts) and current_pts[start_pt_idx][0] < y_min: start_pt_idx += 1 # Accumulate points incrementally as y_max goes higher active_point_ids = [] pt_idx = start_pt_idx for y_max in y_max_choices: while pt_idx < len(current_pts) and current_pts[pt_idx][0] <= y_max: active_point_ids.append(current_pts[pt_idx][1]) pt_idx += 1 # 3. Format output with sequentially generated rectangle IDs rect_id = f"rectangle{rect_counter}" rectangles_output[rect_id] = { "bounds": {"x_min": x_min, "x_max": x_max, "y_min": y_min, "y_max": y_max}, "covered_points": list(active_point_ids) } rect_counter += 1 return rectangles_output, points_dict
GAMS Set Covering Model
$onText
cluster points using axis-aligned rectangles 1. Enumerate all rectangles 2. Set Covering type of model
$offText
*--------------------------------------------------------------- * Settings *---------------------------------------------------------------
option seed=12345;
*--------------------------------------------------------------- * data *---------------------------------------------------------------
scalars maxsize 'size of the outer box' / 10 / k 'number of rectangles' / 5 / ;
sets i 'points' /point1*point50/ r0 'base set rectangles' /rectangle1*rectangle100000/ c 'coordinates' /x,y/ lu 'lower/upper values' /min,max/ ;
parameter p(i,c) 'points'; p(i,c) = uniform(0,maxsize); display p;
*--------------------------------------------------------------- * enumeration algorithm *---------------------------------------------------------------
scalars n 'number of points' nr 'number of possible rectangles' ; n = card(i);
sets r(r0) 'rectangles' contain(i,r0) 'point i is in rect r' ;
Parameter rect(r0,c,lu) 'rectangle coordinates' ;
* create dense data structure for export to Python p(i,c)$(p(i,c)=0) = EPS;
EmbeddedCode Python:
# # number of points #
n = int(list(gams.get("n"))[0])
# # get points and store as list of tuples (id,x,y) #
pgen = gams.get("p",keyFormat=KeyFormat.FLAT) points=[] for i in range(n): px = next(pgen) py = next(pgen) points.append((px[0],px[2],py[2]))
# # import enumeration code #
import sys sys.path.insert(0,'')
import renumerate as renum
print(f"renumerate {renum.version}")
# # start working.... # print(f"points:{n}") rects, mapped_points = renum.find_all_boundary_rectangles(points) print(f"rectangles:{len(rects)}")
# # export to GAMS # gams.set('nr',[len(rects)]) gams.set('r',list(rects.keys())) gams.set('contain', [(p,r) for r in rects.keys() for p in rects[r]['covered_points']] )
coord = [] for r in rects: bounds = rects[r]['bounds'] coord.append((r,'x','min',bounds['x_min'])) coord.append((r,'x','max',bounds['x_max'])) coord.append((r,'y','min',bounds['y_min'])) coord.append((r,'y','max',bounds['y_max'])) gams.set('rect',coord)
endEmbeddedCode nr,r,contain,rect
abort$(nr>card(r0)) "set r0 is too small to store all rectangles";
display nr;
*--------------------------------------------------------------- * calculate area and perimeter *---------------------------------------------------------------
Parameters Area(r0) Perimeter(r0) ;
Area(r) = prod(c,rect(r,c,'max')-rect(r,c,'min')); Perimeter(r) = 2*sum(c,rect(r,c,'max')-rect(r,c,'min'));
*--------------------------------------------------------------- * optimization *---------------------------------------------------------------
binary variable select(r0) 'rectangle selection' ; variables totalArea 'objective' totalPerimeter 'objective' ;
Equations obj1 obj2 cover(i) count ;
obj1.. totalArea =e= sum(r, Area(r)*select(r)); obj2.. totalPerimeter =e= sum(r, Perimeter(r)*select(r));
cover(i).. sum(contain(i,r), select(r)) =g= 1;
count.. sum(r, select(r)) =e= k;
model m1 /all-obj2/; model m2 /all-obj1/;
solve m1 using mip minimizing totalArea; |
K-Medians in R
For comparison, here is an R session that performs a k-medians clustering on our data.
df <- read.table(text=" i x y point1 4.97311132 9.63862756 point2 8.73922492 8.44273500 point3 1.91226253 4.94526579 point4 2.34553216 3.80983451 point5 6.00887463 1.51476321 point6 5.07716895 4.63323882 point7 3.59064347 1.52625715 point8 3.66418827 8.26045475 point9 0.23893508 0.59261597 point10 2.84759145 1.14012489 point11 8.90454986 3.74340300 point12 2.34069434 1.47696552 point13 0.19119162 0.78401936 point14 9.27971966 5.55039788 point15 6.16361826 3.60477909 point16 5.21416943 4.32016162 point17 7.49929678 8.17376279 point18 0.16037081 3.72949918 point19 6.97091674 7.66895295 point20 8.87457833 8.11673997 point21 7.89650246 8.03634493 point22 0.16576575 3.93244606 point23 9.60138480 7.21520825 point24 8.14734411 8.37501191 point25 7.74541956 5.28258128 point26 6.55663865 3.92384422 point27 4.85091115 4.02642582 point28 4.34623796 1.93472329 point29 1.28959803 4.87601902 point30 7.91163148 6.85615242 point31 3.02621818 0.00063517 point32 4.46278421 8.14820541 point33 2.08022534 3.35953314 point34 2.81600217 1.56499252 point35 6.20283636 8.59944445 point36 6.15258941 3.85719583 point37 7.63343999 7.31506452 point38 0.79832230 0.82942221 point39 0.19238302 2.45380498 point40 3.34088869 3.03838803 point41 8.34902953 3.86320590 point42 0.53110352 5.05611011 point43 1.91904761 2.34489699 point44 9.34705030 7.86066922 point45 9.43453787 1.71763581 point46 1.81723665 2.02565095 point47 5.32342079 5.06474382 point48 2.10786534 6.38173925 point49 2.95917598 5.53252693 point50 0.34180259 8.07512765 ",header=T) mat <- as.matrix(df[, c("x", "y")]) library(Kmedians) res <- Kmedians(mat,nclust=5) df$cluster = res$bestresult$cluster library(ggplot2) library(dplyr) # Convert cluster to a factor for discrete coloring df$cluster <- as.factor(df$cluster) # ---------------------------------------------------- # Scatter plot with bounding rectangles # ---------------------------------------------------- # Calculate min/max coordinates for each cluster bounding box cluster_bounds <- df %>% group_by(cluster) %>% summarize( xmin = min(x), xmax = max(x), ymin = min(y), ymax = max(y), .groups = "drop" ) plot <- ggplot() + # Draw the bounding rectangles first so they sit behind the points geom_rect( data = cluster_bounds, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = cluster), alpha = 0.15, # Makes the rectangles semi-transparent color = "black", # Border color linetype = "dashed" ) + # Overlay the points geom_point(data = df, aes(x = x, y = y, color = cluster), size = 3) + labs( title = "Clusters with Bounding Rectangles", x = "X Coordinate", y = "Y Coordinate", color = "Cluster", fill = "Cluster" ) + theme_minimal() print(plot)
The total area is 54.9 and the sum of the perimeters is 68.4.
It is noted that the rectangles are just added for illustration. The underlying mathematical model for the k-medians method is not using them.
No comments:
Post a Comment