👨‍🏭🤖

I just returned from the 2025 INFORMS Computing Society conference, where I had the privilege of organizing a cluster on optimization solvers. The cluster had two sessions, Solvers I and Solvers II, and focussed on new developments in the implementation of optimization solvers. In the coming days, I’m going to explore some of these solvers in more depth. For now, I wanted to give a few hot takes from the sessions while they are still fresh in my mind. ...

In the last post, we used Gurobi’s hierarchical optimization features to compute the Pareto front for primary and secondary objectives in an assignment problem. This relied on Gurobi’s setObjectiveN method and its internal code for managing hierarchical problems. Some practitioners may need to do this without access to a commercial license. This post adapts the previous example to use HiGHS and its native Python interface, highspy. It’s also useful to see what the procedure is in order to understand it better. This isn’t exactly what I’d call hard, but it is easy to mess up.1 ...

One of the first technology choices to make when setting up an optimization stack is which modeling interface to use. Even if we restrict our choices to Python interfaces for MIP modeling, there are lots of options to consider. If you use a specific solver, you can opt for its native Python interface. Examples include libraries like gurobipy, Fusion, highspy, or PySCIPOpt. This approach provides access to important solver-specific features such as lazy constraints, heuristics, and various solver settings. However, it can also lock you into a solver before ready for that. ...

At a Nextmv tech talk a couple weeks ago, I showed a least absolute deviations (LAD) regression model using OR-Tools. This isn’t new – I pulled the formulation from Rob Vanderbei’s “Local Warming” paper, and I’ve shown similar models at conference talks in the past using other modeling APIs and solvers. There are a couple reasons I keep coming back to this problem. One is that it’s a great example of how to build a machine learning model using an optimization solver. Unless you have an optimization background, it’s probably not obvious you can do this. Building a regression or classification model with a solver directly is a great way to understand the model better. And you can customize it in interesting ways, like adding epsilon insensitivity. ...

I attended DPSOLVE 2023 recently and found lots of good inspiration for the next version of Nextmv’s Decision Diagram (DD) solver, Hop. It’s a few years old now, and we learned a lot applying it in the field. Hop formed the basis for our first routing models. While those models moved to a different structure in our latest routing code, the first version broke ground combining DDs with Adaptive Large Neighborhood Search (ALNS), and its use continues to grow organically. ...

Note: This post has been updated to work with HiGHS. A fun geometry problem to think about is: given two polygons, do they intersect? That is, do they touch on the border or overlap? Does one reside entirely within the other? While this question has obvious applications in computer graphics (see: arcade games of the 1980s), it’s also important in areas such as cutting and packing problems. There are a number of way to answer this. In computer graphics, the problem is often approached using a clipping algorithm. This post examines a couple of simpler techniques using linear inequalities and properties of convexity. To simplify the presentation, we assume we’re only interested in convex polygons in two dimensions. We also assume that rotation is not an issue. That is, if one of the polygons is rotated, we can simply re-test to see if they overlap. ...

Note: This post was originally written using Julia v0.2, GLPK, and Hedonometer data through 2014. It has been updated to use Julia v1.11, HiGHS, and data through May 26, 2025. Hedonometer popped onto my radar a couple weeks ago. It’s a nifty project, attempting to convert samples of words found in the Twitter Gardenhose feed into a time series of happiness. While I’m not a computational social scientist, I must say the data does have a nice intuitive quality to it. There are obvious trends in happiness associated with major holidays, days of the week, and seasons. It seems like the sort of data that could be decomposed into trends based on those various components. The Hedonometer group has, of course, done extensive analyses of their own data which you can find on their papers page. ...

In the previous post, we considered preprocessing for the vehicle routing problem where the vehicles have different starting locations. Our goal was to create potentially overlapping regions for the entire US which we could later use for route construction. We defined these regions using all 5-digit zip codes in the continental US for which one of our regional headquarters is the closest, or one of $n$ closest, headquarters in terms of Euclidean distance. The resulting regions gave us some flexibility in terms of how much redundancy we allow in our coverage of the country. ...

Consider an instance of the vehicle routing problem in which we have drivers that are geographically distributed, each in a unique location. Our goal is to deliver goods or services to a set of destinations at the lowest cost. It does not matter to our customers which driver goes to which destination, so long as the deliveries are made. One can think of this problem as a collection of travelling salesman problems, where there are multiple salespeople in different locations and a shared set of destinations. We attempt to find the minimum cost schedule for all salespeople that visits all destinations, where each salesman can technically go anywhere. ...

Note: This post was written before Gurobi supported nonlinear optimization. It has been updated to work with Python 3. A common problem in handling geometric data is determining the center of a given polygon. This is not quite so easy as it sounds as there is not a single definition of center that makes sense in all cases. For instance, sometimes computing the center of a polygon’s bounding box may be sufficient. In some instances this may give a point on an edge (consider a right triangle). If the given polygon is non-convex, that point may not even be inside or on its boundary. ...

Note: A reader pointed out that Union-Find is a very efficient way to accomplish this task. Start there if you have the same problem! Last week, Paul Rubin wrote an excellent post on Extracting a Connected Graph from an existing graph. Lately I’ve been performing related functions on data from OpenStreetMap, though without access to a solver. In my case I’m taking in arbitrary network data and splitting it into disconnected sub-networks. I thought it might be a good case study to show an algorithmic way doing this and some of the performance issues I ran into. ...

Note: This post was updated to work with Python 3 and the 2nd edition of “Integer Programming” by Laurence Wolsey. We’ve been studying Lagrangian Relaxation (LR) in the Advanced Topics in Combinatorial Optimization course I’m taking this term, and I had some difficulty finding a simple example covering its application. In case anyone else finds it useful, I’m posting a Python version for solving the Generalized Assignment Problem (GAP). This won’t discuss the theory of LR at all, just give example code using Gurobi. ...

Note: This post was updated to work with Python 3 and PySCIPOpt. The original version used Python 2 and python-zibopt. It has also been edited for clarity. As a followup to the last post, I created another SCIP example for finding Normal Magic Squares. This is similar to solving a Sudoku problem, except that here the number of binary variables depends on the square size. In the case of Sudoku, each cell has 9 binary variables – one for each potential value it might take. For a normal magic square, there are $n^2$ possible values for each cell, $n^2$ cells, and one variable representing the row, column, and diagonal sums. This makes a total of $n^4$ binary variables and one continuous variables in the model. ...

Note: This post was updated to work with Python 3 and PySCIPOpt. The original version used Python 2 and python-zibopt. It has also been edited for clarity. Back in October of 2011, I started toying with a model for finding magic squares using SCIP. This is a fun modeling exercise and a challenging problem. First one constructs a square matrix of integer-valued variables. from pyscipopt import Model # [...snip...] m = Model() matrix = [] for i in range(size): row = [m.addVar(vtype="I", lb=1) for _ in range(size)] for x in row: m.addCons(x <= M) matrix.append(row) Then one adds the following constraints: ...

In response to this post, Ben Bitdiddle inquires: I understand the concept of using a companion set to remove duplicates from a list while preserving the order of its elements. But what should I do if these elements are composed of smaller pieces? For instance, say I am generating combinations of numbers in which order is unimportant. How do I make a set recognize that [1,2,3] is the same as [3,2,1] in this case? ...

This is based on a lightning talk I gave at the LA PyLadies October Hackathon. I’m actually not going to go into anything much resembling algorithmic complexity here. What I’d like to do is present a common performance anti-pattern that I see from novice programmers about once every year or so. If I can prevent one person from committing this error, this post will have achieved its goal. I’d also like to show how an intuitive understanding of time required by operations in relation to the size of data they operate on can be helpful. ...

I find I have to build simulations with increasing frequency in my work and life. Usually this indicates I’m faced with one of the following situations: The need for a quick estimate regarding the quantitative behavior of some situation. The desire to verify the result of a computation or assumption. A situation which is too complex or random to effectively model or understand. Anyone familiar at all with simulation will recognize the last item as the motivating force of the entire field. Simulation models tend to take over when systems become so complex that understanding them is prohibitive in cost and time or entirely infeasible. In a simulation, the modeler can focus on individual interactions between entities while still hoping for useful output in the form of descriptive statistics. ...

Note: This post was updated to work with NetworkX and for clarity. It’s possible this will turn out like the day when Python 2.5 introduced coroutines. At the time I was very excited. I spent several hours trying to convince my coworkers we should immediately abandon all our existing Java infrastructure and port it to finite state machines implemented using Python coroutines. After a day of hand waving over a proof of concept, we put that idea aside and went about our lives. ...

I recently stumbled across an implementation of the affine scaling interior point method for solving linear programs that I’d coded up in R once upon a time. I’m posting it here in case anyone else finds it useful. There’s not a whole lot of thought given to efficiency or numerical stability, just a demonstration of the basic algorithm. Still, sometimes that’s exactly what one wants. solve.affine <- function(A, rc, x, tolerance=10^-7, R=0.999) { # Affine scaling method while (T) { X_diag <- diag(x) # Compute (A * X_diag^2 * A^t)-1 using Cholesky factorization. # This is responsible for scaling the original problem matrix. q <- A %*% X_diag**2 %*% t(A) q_inv <- chol2inv(chol(q)) # lambda = q * A * X_diag^2 * c lambda <- q_inv %*% A %*% X_diag^2 %*% rc # c - A^t * lambda is used repeatedly foo <- rc - t(A) %*% lambda # We converge as s goes to zero s <- sqrt(sum((X_diag %*% foo)^2)) # Compute new x x <- (x + R * X_diag^2 %*% foo / s)[,] # If s is within our tolerance, stop. if (abs(s) < tolerance) break } x } This function accepts a matrix A which contains all technological coefficients for an LP, a vector rc containing its reduced costs, and an initial point x interior to the LP’s feasible region. Optional arguments to the function include a tolerance, for detecting when the method is within an acceptable distance from the optimal point, and a value for R, which must be strictly between 0 and 1 and controls scaling. ...

Note: This post was edited for clarity. For the final JAPH in this series, I implemented a simple transpiler that converts a small subset of Scheme programs to equivalent Python programs. It starts with a Scheme program that prints 'just another scheme hacker'. (define (output x) (if (null? x) "" (begin (display (car x)) (if (null? (cdr x)) (display "\n") (begin (display " ") (output (cdr x))))))) (output (list "just" "another" "scheme" "hacker")) The program then tokenizes that Scheme source, parses the token stream, and converts that into Python 3. ...

Note: This post was edited for clarity. This JAPH uses a Turing machine. The machine accepts any string that ends in '\n' and allows side effects. This lets us print the value of the tape as it encounters each character. While the idea of using lambda functions as side effects in a Turing machine is a little bizarre on many levels, we work with what we have. And Python is multi-paradigmatic, so what the heck. ...

Note: This post was edited for clarity. At this point, tricking python into printing strings via indirect means got a little boring. So I switched to obfuscating fundamental computer science algorithms. Here’s a JAPH that takes in a Huffman coded version of 'just another python hacker', decodes, and prints it. # Build coding tree def build_tree(scheme): if scheme.startswith('*'): left, scheme = build_tree(scheme[1:]) right, scheme = build_tree(scheme) return (left, right), scheme else: return scheme[0], scheme[1:] def decode(tree, encoded): ret = '' node = tree for direction in encoded: if direction == '0': node = node[0] else: node = node[1] if isinstance(node, str): ret += node node = tree return ret tree = build_tree('*****ju*sp*er***yct* h**ka*no')[0] print( decode(tree, bin(10627344201836243859174935587).lstrip('0b').zfill(103)) ) The decoding tree is like a LISP-style sequence of pairs. '*' represents a branch in the tree while other characters are leaf nodes. This looks like the following. ...

Note: This post was updated to work with Python 3.12. It may not work with different versions. Here’s a JAPH composed solely for effect. For each letter in 'just another python hacker' it loops over each the characters ' abcdefghijklmnopqrstuvwxyz', printing each. Between characters it pauses for 0.05 seconds, backing up and moving on to the next if it hasn’t reached the desired one yet. This achieves a sort of rolling effect by which the final string appears on our screen over time. ...

Note: This post was updated to work with Python 3.12. It may not work with different versions. No series of JAPHs would be complete without ROT13. This is the example through which aspiring Perl programmers learn to use tr and its synonym y. In Perl the basic ROT13 JAPH starts as: $foo = 'whfg nabgure crey unpxre'; $foo =~ y/a-z/n-za-m/; print $foo; Python has nothing quite so elegant in its default namespace. However, this does give us the opportunity to explore a little used aspect of strings: the translate method. If we construct a dictionary of ordinals we can accomplish the same thing with a touch more effort. ...

Here’s the second in my reformed JAPH series. It takes an anagram of 'just another python hacker' and converts it prior to printing. It sorts the anagram by the indices of another string, in order of their associated characters. This is sort of like a pre-digested Schwartzian transform. x = 'upjohn tehran hectors katy' y = '1D0HG6JFO9P5ICKAM87B24NL3E' print(''.join(x[i] for i in sorted(range(len(x)), key=lambda p: y[p]))) Obfuscation consists mostly of using silly machinations to construct the string we use to sort the anagram. ...

Note: This post was edited for clarity. Many years ago, I was a Perl programmer. Then one day I became disillusioned at the progress of Perl 6 and decided to import this. This seems to be a fairly common story for Perl to Python converts. While I haven’t looked back much, there are a number of things I really miss about perl (lower case intentional). I miss having value types in a dynamic language, magical and ill-advised use of cryptocontext, and sometimes even pseudohashes because they were inexcusably weird. A language that supports so many ideas out of the box enables an extended learning curve that lasts for many years. “Perl itself is the game.” ...

I hope you saw “China’s way to the top” on the Post’s website recently. It’s a very clear presentation of their statement and is certainly worth a look. So say you’re an economist and you actually do need to produce a realistic estimate of when China’s GDP surpasses that of the USA. Can you use such an approach? Not really. There are several simplifying assumptions the Post made that are perfectly reasonable. However, if the goal is an analytical output from a highly random system such as GDP growth, one should not assume the inputs are fixed. (I’m not saying I have any gripe with their interactive. This post has a different purpose.) ...

Note: This post was updated to include an example data file. I thought it might be useful to follow up the last post with another one showing the same examples in R. R provides a function called lm, which is similar in spirit to NumPy’s linalg.lstsq. As you’ll see, lm’s interface is a bit more tuned to the concepts of modeling. We begin by reading in the example CSV into a data frame: ...

This post is based on a memo I sent to some former colleagues at the Post. I’ve edited it for use here since it fits well as the second in a series on simple data fitting techniques. If you’re among the many enlightened individuals already using regression analysis, then this post is probably not for you. If you aren’t, then hopefully this provides everything you need to develop rudimentary predictive models that yield surprising levels of accuracy. ...

Consider this scenario: You run a contest that accepts votes from the general Internet population. In order to encourage user engagement, you record any and all votes into a database over several days, storing nothing more than the competitor voted for, when each vote is cast, and a cookie set on the voter’s computer along with their apparent IP addresses. If a voter already has a recorded cookie set they are denied subsequent votes. This way you can avoid requiring site registration, a huge turnoff for your users. Simple enough. ...

Note: This post was updated to work with Python 3 and PySCIPOpt. The original version used Python 2 and python-zibopt. Data fitting is one of those tasks that everyone should have at least some exposure to. Certainly developers and analysts will benefit from a working knowledge of its fundamentals and their implementations. However, in my own reading I’ve found it difficult to locate good examples that are simple enough to pick up quickly and come with accompanying source code. ...

Note: This post was updated to work with Python 3. One of the most useful tools one learns in an Operations Research curriculum is Monte Carlo Simulation. Its utility lies in its simplicity: one can learn vital information about nearly any process, be it deterministic or stochastic, without wading through the grunt work of finding an analytical solution. It can be used for off-the-cuff estimates or as a proper scientific tool. All one needs to know is how to simulate a given process and its appropriate probability distributions and parameters if that process is stochastic. ...

One of the difficulties we encounter in solving the Traveling Salesman Problem (TSP) is that, for even a small number of cities, a complete description of the problem requires a factorial number of constraints. This is apparent in the standard formulation used to teach the TSP to OR students. Consider a set of $n$ cities with the distance from city $i$ to city $j$ denoted $d_{ij}$. We attempt to minimize the total distance of a tour entering and leaving each city exactly once. $x_{ij} = 1$ if the edge from city $i$ to city $j$ is included in the tour, $0$ otherwise: ...

Uncapacitated Lot Sizing (ULS) is a classic OR problem that seeks to minimize the cost of satisfying known demand for a product over time. Demand is subject to varying costs for production, set-up, and storage of the product. Technically, it is a mixed binary integer linear program – the key point separating it from the world of linear optimization being that production cannot occur during any period without paying that period’s fixed costs for set-up. Thus it has linear nonnegative variables for production and storage amounts during each period, and a binary variable for each period that determines whether or not production can actually occur. ...