The Power of Constraint Solvers

A developer with a career in constraint solving and constrained optimization surveys classical AI techniques — SAT, SMT, constraint programming, and derivative-free/surrogate-based optimization — and demonstrates them by building a minimal-depth 1-bit full-adder circuit from NAND gates. The writeup walks through real-world applications including verifying distributed program specifications for race conditions and deadlocks, generating mechatronic system architectures, and tuning bicycle aerodynamics and boiler modulation ratios.

With my kids 1 footnote-1 having derailed every single one of my side projects including this blog , I thought to revisit a topic very dear to me: constraint solving and constrained optimization . I already wrote an article about “ Machine Reasoning https://btmc.substack.com/p/machine-reasoning-the-forgotten-side ” 2 footnote-2 , which is a made up umbrella term for all sorts of classical AI techniques in contrast to modern Machine Learning https://en.wikipedia.org/wiki/Machine learning -based approaches , but it was a very high level overview focused on listing different types of constraint solvers that didn’t really provide any particularly useful piece of information. Not very interesting, so here’s a better one. My entire professional career and academic research has involved constraint solving or constrained optimization in some capacity. Here are some use cases the tools I worked on have been applied to: - Verifying that distributed program specifications were free of race conditions and deadlocks, and then generating correct-by-construction program skeletons containing just the network communication from those specifications. - Generating mechatronic system architectures, including electrical power systems in aircraft and gearboxes in cars. - Parameter optimization, including improving bicycle aerodynamics and greatly improving the modulation ratio https://homesteadenergy.co.uk/what-is-boiler-modulation-and-how-does-it-save-money/ of a domestic water boiler. Each of the three bullet points used a different type of solver. The first used a chain of verification tools that ultimately invoked a Satisfiability Modulo Theories https://en.wikipedia.org/wiki/Satisfiability modulo theories SMT solver. The second used a combination of a Boolean Satisfiability https://en.wikipedia.org/wiki/Boolean satisfiability problem SAT solver and a Constraint Programming https://en.wikipedia.org/wiki/Constraint programming CP solver. The last one used a custom solver implementing a large portfolio of Derivative Free https://en.wikipedia.org/wiki/Derivative-free optimization and Surrogate-based https://en.wikipedia.org/wiki/Surrogate model optimization algorithms. Most programmers are completely unfamiliar with these sorts of tools and the wide variety of very hard problems they can tackle. If you are into game dev you might be familiar with Wave Function Collapse https://github.com/mxgmn/WaveFunctionCollapse , which is basically a Constraint Programming technique https://en.wikipedia.org/wiki/AC-3 algorithm applied to the procedural generation of textures and tile-maps. All those LLM math proofs coming out are using Lean https://lean-lang.org a proof assistant and are probably making very liberal use of the ‘ Grind https://lean-lang.org/doc/reference/latest/The--grind--tactic/ ’ 3 footnote-3 tactic for automated proof search, which uses SMT solver techniques under the hood. But this is all a little too abstract, so let’s solve a cool problem together. The Problem: An Adder Circuit Here’s what we want to accomplish: Using nothing but NAND-gates https://en.wikipedia.org/wiki/NAND gate , a “functionally complete” logic gate that can be used to implement any boolean formula, implement a 1-bit Full-Adder <https://en.wikipedia.org/wiki/Adder electronics Full adder circuit. A NAND-gate takes 2 input bits A and B and outputs one bit C that is the negation of the conjunction of A and B: A Full-Adder takes 2 input bits A and B plus an input carry bit C<sub i</sub , and outputs a sum bit S plus an output carry bit C<sub o</sub 4 footnote-4 , behaving as follows: So for example, if A=1, B=0, and C<sub i</sub =1, then the output is C<sub o</sub =1 and S=0. Your CPU and GPU are full of adder circuits pun intended , though not made up of only NAND-gates. We want to use the bare minimum amount of NAND-gates necessary, and we want the configuration with the least depth as an approximation of propagation delay less delay = better performance . Here’s a solution: 9 NAND gates. Depth 6. Would you be able to come up with this solution by hand? Very likely, it’s not a hard problem. But can you write a program that finds it ? Also, is this the best we can do? If we increase the depth can we get away with less gates? If we use more gates can we pull off a lower depth? Or maybe can we improve both metrics at the same time? Give that some thought before proceeding with the rest of the article. Constraint Satisfaction Problems The problem above, of finding that Full-Adder circuit, can be modeled as a Constraint Satisfaction Problem https://en.wikipedia.org/wiki/Constraint satisfaction problem CSP . Any NP problem <https://en.wikipedia.org/wiki/NP complexity can be modeled as a CSP because constraint satisfaction is itself an NP-Complete https://en.wikipedia.org/wiki/NP-completeness problem. If we add the objective of finding the lowest amount of NAND-gates necessary and/or the lowest depth required, then it becomes a Constrained Optimization Problem https://en.wikipedia.org/wiki/Constrained optimization COP . Multi-objective optimization is a can of worms 5 footnote-5 , so we’ll do a little trick to have a single objective: For a given fixed number of NAND-gates, find the lowest depth. If there is no solution, increase the number of NAND-gates. Unless P=NP https://en.wikipedia.org/wiki/P versus NP problem , any algorithm that can solve this problem must have worst-case exponential complexity: O 2<sup n</sup . That might sound really bad, but it’s worst-case complexity. In practice, state-of-the-art solvers can eat a problem like the above for breakfast, as you’ll soon see. Modeling the Problem A constraint satisfaction problem consists of decision variables, the domains of those variables what values they can take , constraints relating those variables, and in the case of constrained optimization, an objective. We’ll make use of the excellent Python package CPMpy https://cpmpy.readthedocs.io/en/latest/ because it has a super nice API and can invoke a ton of different solvers under the hood. Here’s our starting point: python import cpmpy as cp def main : model = cp.Model We'll model the problem here solver = cp.SolverLookup.get "ortools", model has solution = solver.solve print f"Status: {solver.status }" if has solution: We'll output the solution here pass if name == " main ": main OR-Tools https://developers.google.com/optimization is Google’s suite of constraint solving and optimization software. The solver lookup above is actually loading a specific solver from that suite, CP-SAT https://developers.google.com/optimization/cp/cp solver , which has won a constraint solving competition https://www.minizinc.org/challenge/ every year for 13 years straight, it’s that good. If you run this program you’ll get the following result: Status: ExitStatus.FEASIBLE 0.019812 seconds There are no variables and no constraints, so the problem is trivial. Other statuses we can get include: - OPTIMAL: The solver found a solution and proved it is the best possible one. - UNSATISFIABLE: The solver proved there is no solution. - UNKNOWN: The solver was not able to solve the problem in the time given. Now let us start setting up the problem. First, we’ll just model the inputs and outputs of the system and impose that it does indeed add the bits: a = cp.boolvar name="a" b = cp.boolvar name="b" c i = cp.boolvar name="c i" s = cp.boolvar name="s" c o = cp.boolvar name="c o" model.add 2 c o + s == a + b + c i CPMpy uses a lot of operator overloading abuse to let you write the above If we run the program we’ll get FEASIBLE again and it’ll take about the same time. We can also print the solution it found within if has solution : if has solution: print f"{a.name}: {a.value }" print f"{a.name}: {b.value }" print f"{c i.name}: {c i.value }" print f"{s.name}: {s.value }" print f"{c o.name}: {c o.value }" I just get false on all of them, which makes sense, 2 0 + 0 = 0 + 0 + 0. Try adding a constraint setting c o == True and see what happens make sure to remove it after . But we’re just directly relating the outputs to the inputs here, this tells us nothing useful. What we want is to propagate those inputs through some number of NAND gates, such that, for all possible inputs , the output is as expected. We’ll do this in multiple steps. Modeling NAND gates Each NAND gate will be modeled as 4 variables: 2 inputs, 1 output, and one variable tracking the depth where this NAND gate sits: python NUM GATES = 9 class NANDGate: def init self, id: int : self.id = id self.a = cp.boolvar name=f"nand {id} a" self.b = cp.boolvar name=f"nand {id} b" self.c = cp.boolvar name=f"nand {id} c" we need to give lower and upper bounds to integer variables depth can never exceed the number of gates, so max is NUM GATES self.depth = cp.intvar 1, NUM GATES, name=f"nand {id} depth" self.connections to a = self.connections to b = self.connections from c = def constrain value computation self, model : model.add self.c == ~ self.a & self.b The class itself is just a way to bundle the data together, it serves no other purpose. We’ll tackle those connection lists later more variables . The depth will also come into play later. For now, we just want to add some NAND gates to the problem and enforce that they compute the right value, meaning output C has to be the negation of the conjunction of inputs A and B. nand gates = for i in range 0, NUM GATES : gate = NANDGate i+1 nand gates.append gate gate.constrain value computation model We can also print the values of these new variables in the solution, as a sanity check. It’s not really the information we ultimately care about, but it tells us if we screwed something up in the modeling, you can remove this later if you want: if has solution: ... for gate in nand gates: print f"{gate.a.name}: {gate.a.value }" print f"{gate.b.name}: {gate.b.value }" print f"{gate.c.name}: {gate.c.value }" Modeling Connections Now we want these NAND gates to connect to each other. Meaning, if we decide that NAND gate 1’s output C is connected to NAND gate 2’s input A, then the values of C and A must be equal. We’ll model connections as follows: python class GateConnections: def init self, from gate: NANDGate, to gate: NANDGate : self.c a = cp.boolvar name=f"nand {from gate.id} c to nand {to gate.id} a" self.c b = cp.boolvar name=f"nand {from gate.id} c to nand {to gate.id} b" self.from gate = from gate self.to gate = to gate def constrain value propagation self, model : model.add self.c a.implies self.from gate.c == self.to gate.a model.add self.c b.implies self.from gate.c == self.to gate.b def compute connections gates: list NANDGate : connections = for i in range len gates : for j in range i + 1, len gates : con = GateConnections gates i , gates j connections.append con gates j .connections to a.append con.c a gates j .connections to b.append con.c b gates i .connections from c.append con.c a gates i .connections from c.append con.c b return connections Each GateConnections instance is actually modeling 2 different connections between 2 gates, one for each of the input ports of the target. We could have kept them separate but this shortens the blog post. Each of the connections from output to input is modeled as a boolean variable signaling that decision. If c a is true, then that means c is connected to a . Note that an output port can connect to more than one input port, so having both c a and c b be true is perfectly fine. The constrain value propagation method adds constraints that ensure that the value of the output is propagated to the target input if the connection exists. There are no backward connections in an adder, so the loop that creates connections and the respective connection variables only creates them from lower ID gates to higher ID gates. This cuts a lot of the search space that consists of nothing but swaps of gate IDs, a form of symmetry breaking https://en.wikipedia.org/wiki/Symmetry-breaking constraints . Now, the connection constraints: connections = compute connections nand gates for con in connections: con.constrain value propagation model for gate in nand gates: model.add cp.sum gate.connections to a == 1 model.add cp.sum gate.connections to b == 1 model.add cp.sum gate.connections from c = 1 The first loop should be obvious, but for the second loop we’re saying the following: for each NAND gate, each of its input ports must be connected to exactly 1 other port exactly one of the connection variables involving it must be true , and the output port must connect at least once, but can connect any number of times at least one of the connection variables involving it must be true . If you actually try to run the program now, you’ll get ExitStatus.UNSATISFIABLE . The reason is that there’s no way to connect all the ports with 9 NAND gates that can only connect forward: the first NAND gate’s inputs can’t connect to anything and the last NAND gate’s output can’t connect to anything Modeling Mappings We “forgot” pedagogically about the input and output ports of the adder itself. We need to implement a mapping from them to the NAND gates’ ports. We model it similarly to connections, but a mapping is input to input or output to output. First a small refactoring: python class Adder: def init self : self.a = cp.boolvar name="a" self.b = cp.boolvar name="b" self.c i = cp.boolvar name="c i" self.s = cp.boolvar name="s" self.c o = cp.boolvar name="c o" self.a mappings = self.b mappings = self.c i mappings = self.s mappings = self.c o mappings = def constrain value computation self, model : model.add 2 self.c o + self.s == self.a + self.b + self.c i The loose variables that represented the adder are now also bundled together, and we also track their mappings much like we tracked connections in the NANDGate instances. There is only one Adder instance but now we can pass it around. Our main function now starts as follows: model = cp.Model adder = Adder adder.constrain value computation model nand gates = ... With this small refactoring in place, we can implement mappings: python class GateMappings: def init self, gate: NANDGate, adder: Adder : self.gate = gate self.adder = adder self.a a = cp.boolvar name=f"adder a nand {gate.id} a" self.a b = cp.boolvar name=f"adder a nand {gate.id} b" self.b b = cp.boolvar name=f"adder b nand {gate.id} b" self.b a = cp.boolvar name=f"adder b nand {gate.id} a" self.c i a = cp.boolvar name=f"adder c i nand {gate.id} a" self.c i b = cp.boolvar name=f"adder c i nand {gate.id} b" self.s c = cp.boolvar name=f"adder s nand {gate.id} c" self.c o c = cp.boolvar name=f"adder c o nand {gate.id} c" def constrain value propagation self, model : model.add self.a a.implies self.gate.a == self.adder.a model.add self.b b.implies self.gate.b == self.adder.b model.add self.a b.implies self.gate.b == self.adder.a model.add self.b a.implies self.gate.a == self.adder.b model.add self.c i a.implies self.gate.a == self.adder.c i model.add self.c i b.implies self.gate.b == self.adder.c i model.add self.s c.implies self.gate.c == self.adder.s model.add self.c o c.implies self.gate.c == self.adder.c o These mappings act much like connections, and we can treat them as such: python def compute mappings adder: Adder, gates: list NANDGate : mappings = for gate in gates: mapping = GateMappings gate, adder mappings.append mapping gate.connections to a.append mapping.a a gate.connections to a.append mapping.b a gate.connections to b.append mapping.b b gate.connections to b.append mapping.a b gate.connections to a.append mapping.c i a gate.connections to b.append mapping.c i b gate.connections from c.append mapping.s c gate.connections from c.append mapping.c o c adder.a mappings.append mapping.a a adder.a mappings.append mapping.a b adder.b mappings.append mapping.b b adder.b mappings.append mapping.b a adder.c i mappings.append mapping.c i a adder.c i mappings.append mapping.c i b adder.s mappings.append mapping.s c adder.c o mappings.append mapping.c o c return mappings Note that we’re also adding mappings to the lists of connections in the gates. That’s because if an input port is mapped, it cannot also be connected to another port, the constraint that is it connected exactly once must also consider the option of mapping it instead similar logic applies to output ports . Then we can add the constraints: add this before the connection sum constraints so the mappings are included there mappings = compute mappings nand gates, adder for mapping in mappings: mapping.constrain value propagation model The problem is solvable now but the configuration I got most definitely does not compute an addition: What’s going on here? For All Possible Inputs The issue is this constraint: model.add 2 self.c o + self.s == self.a + self.b + self.c i This constraint just says the solver needs to find a value for these variables that makes the constraint true, it doesn’t say that for all possible inputs the correct sum is computed. So it can simply pick 1 input/output pair where the configuration of NAND gates happens to give the right result, even if the result is wrong for all other inputs. Imagine we were doing something much simpler: implementing AND. When both inputs are True or both inputs are False, AND and OR agree, so an implementation of “OR” would pass as a solution for “AND” for certain inputs. Does this make sense? We need it to give the right solution for all inputs, simultaneously. That “for all” there in the sentence above is quite important, it is the universal quantifier https://en.wikipedia.org/wiki/Universal quantification ∀ . The problem we’re trying to solve is actually of the form: Which is a synthesis problem, it states: there exists some program p such that, for all inputs i , the program constraints P p imply the specification constraints S p, i for those inputs. There are various ways to deal with that universal quantifier: - Some solvers have native support for quantifiers, but performance tends to be rather dubious it’s a very hard problem in general . - Two solvers can work together in a Counter-Example Guided Inductive Synthesis https://en.wikipedia.org/wiki/Program synthesis Counter-example guided inductive synthesis CEGIS loop, each handling one of the quantifiers. - We can just unroll the universal quantifier for all possible inputs. We only have 8 possible inputs, so the latter option is easily the best choice for this Full-Adder use case. But how do we actually model that? Each variable can only take one value, and we want to propagate 8 different inputs through them. The answer is to create 8 different sets of the “value” variables, and each input will be routed through one of those sets. The connection variables will then route all 8 sets of value variables simultaneously. We’ll start by making the following changes: python class NANDGate: def init self, id: int : self.id = id Note the 8 as the first input, a is now an array of 8 boolean variables self.a = cp.boolvar 8, name=f"nand {id} a" self.b = cp.boolvar 8, name=f"nand {id} b" self.c = cp.boolvar 8, name=f"nand {id} c" ... class Adder: def init self : self.a = cp.boolvar 8, name="a" self.b = cp.boolvar 8, name="b" self.c i = cp.boolvar 8, name="c i" self.s = cp.boolvar 8, name="s" self.c o = cp.boolvar 8, name="c o" All of the value variables get an 8 as the first argument, turning them into arrays. CPMpy has an array-oriented API, so a + b is the same as a 0 + b 0 , a 1 + b 1 , a 2 + b 2 , … . We do need to change the connection implication constraints because they don’t work directly on arrays, we need to say that we want the conjunction of those expressions, using cp.all : python class GateConnections: ... def constrain value propagation self, model : model.add self.c a.implies cp.all self.from gate.c == self.to gate.a model.add self.c b.implies cp.all self.from gate.c == self.to gate.b class GateMappings: ... def constrain value propagation self, model : model.add self.a a.implies cp.all self.gate.a == self.adder.a model.add self.b b.implies cp.all self.gate.b == self.adder.b model.add self.a b.implies cp.all self.gate.b == self.adder.a model.add self.b a.implies cp.all self.gate.a == self.adder.b model.add self.c i a.implies cp.all self.gate.a == self.adder.c i model.add self.c i b.implies cp.all self.gate.b == self.adder.c i model.add self.s c.implies cp.all self.gate.c == self.adder.s model.add self.c o c.implies cp.all self.gate.c == self.adder.c o Now, we’re ready for the input constraints, we just need to define the truth table: python def constrain truth table model, adder: Adder : table = 0, 0, 0, 0, 0 , 1, 0, 0, 1, 0 , 0, 1, 0, 1, 0 , 0, 0, 1, 1, 0 , 1, 1, 0, 0, 1 , 1, 0, 1, 0, 1 , 0, 1, 1, 0, 1 , 1, 1, 1, 1, 1 , for i in range 8 : array = cp.cpm array adder.a i , adder.b i , adder.c i i , adder.s i , adder.c o i model.add array == table i def main : ... right before we lookup the solver: constrain truth table model, adder We now have a working full-adder generator. We didn’t write a visualizer so it’s hard to tell on your end I have a vibe coded one that exports the solution to netlistsvg https://github.com/nturley/netlistsvg , but we could also write an exporter to some logic circuit emulator and test it there. Anyway, just believe me for now. Run the program and it’ll very quickly tell you the problem is feasible and provide a solution. Here’s what mine looks like: Looks rather… familiar? It’s the same as the one we were looking for, just shifted around. Running the program a bunch of times seems to generate very similar images each time, but that’s no confirmation. If we ask for all solutions using solver.solveAll instead of solver.solve , we get that there are 24576 solutions, but many of those solutions are symmetric meaning equivalent , they just have elements swapped around that don’t make any difference to the result. Dealing with symmetries and graph isomorphisms https://en.wikipedia.org/wiki/Graph isomorphism is out of scope for this article, however. But here’s the million dollar question: Is there a way to make a Full-Adder with only 8 NAND-gates? We just have to set NUM GATES to 8 and solve again. In just 99 milliseconds on my dying M1 MacBook Pro I’m told in no uncertain terms that nope, you cannot make a Full-Adder with only 8 NAND-gates: Status: ExitStatus.UNSATISFIABLE 0.099592 seconds As a reminder, this doesn’t mean the solver failed to find a solution, it explored the entire design space in 99 milliseconds. There is no solution . Some solvers can even produce machine-verifiable albeit unreadable mathematical proofs that there really is no solution. If we remove the truth table constraints we get nonsensical solutions again, but now with 8 gates instead of 9. What about depth though? Can we do better than 6? Modeling Depth Spoiler alert: no, 6 is the best we can do, there is in fact only 1 way to make a Full-Adder out of NAND-gates those 24576 solutions are all equivalent . But let’s get confirmation anyway. We already added the depth variables, but we didn’t relate them in any way. The easiest way to fit depth into our existing constraints is to enforce that if two gates are connected, then the depth of the target gate is greater: python def constrain depth model, connections: list GateConnections : pairs = set for con in connections: if con.from gate.id, con.to gate.id in pairs: continue model.add con.c a.implies con.to gate.depth con.from gate.depth model.add con.c b.implies con.to gate.depth con.from gate.depth pairs.add con.from gate.id, con.to gate.id Then our objective is to minimize the max depth, simple enough: python def main : ... constrain depth model, connections model.minimize cp.max gate.depth for gate in nand gates ... if has solution: print f"Depth: {model.objective value }" And the result is, sadly: Status: ExitStatus.OPTIMAL 0.476214 seconds Depth: 6 Depth 6 is as low as we can go with NAND gates. Neat How do I implement a Solver? This is normally the part of the blog post where I tell you how to implement some crazy thing in C, in this case a constraint solver, but sadly this time I won’t. The article is already huge and a competitive constraint solver requires quite a lot of engineering. It wouldn’t be an article, more like a book. Here are some starting points for different solver types: - Conflict-Driven Clause Learning https://en.wikipedia.org/wiki/Conflict-driven clause learning CDCL, used in SAT . - Arc-Consistency Algorithm 3 https://en.wikipedia.org/wiki/AC-3 algorithm AC-3, used in CP - Simplex Algorithm https://en.wikipedia.org/wiki/Simplex algorithm used in Linear Programming https://en.wikipedia.org/wiki/Linear programming . Maybe I’ll make a “Crafting Solvers” in honor of Bob Nystrom’s “ Crafting Interpreters https://craftinginterpreters.com ” someday. Conclusion Constraint Solving is an extremely powerful approach for tackling hard problems that would otherwise require highly-specialized algorithms and heuristics. Despite having exponential complexity in general, state-of-the-art solvers can chew through many real-world problems with ease. Sadly with larger adders it begins to struggle 2-bit is still fair game , though you can always make larger adders out of Full-Adders as opaque blocks to make it tractable though the optimal solution then involves adding extra gates to handle the carry computations in parallel, to minimize propagation delay . Here’s a gist https://gist.github.com/sirwhinesalot/a05f772e1d0b8a39902af90b8a69c890 with the complete code above. See you on the next one. 1 footnote-anchor-1 Totally not my atrocious time management skills, it is 100% the kids’ fault. 2 footnote-anchor-2 Not to be confused with “reasoning” in LLMs. 3 footnote-anchor-3 The Isabelle theorem prover calls the same sort of tactic ‘sledgehammer’. Now that’s a name. 4 footnote-anchor-4 A half-adder would lack the input carry. 5 footnote-anchor-5 Many solvers do not natively support multiple objectives, and there are many different ways of setting up the problem when there are multiple competing objectives. Further reading: Multi-objective Optimization https://en.wikipedia.org/wiki/Multi-objective optimization Solution