hey saw your question about a solver for systems of non-linear inequalities classic problem this one and it's a rabbit hole for sure the first thing you have to get is that this is not like solving `Ax = b` this isn't high school algebra where you find a nice single point or a clean line where things intersect with non-linear stuff you're often looking for a *region* a feasible region and that region can be disconnected it can be a weird blob it can have holes in it it can be non-convex which is a huge pain and finding just a single point inside that region can be computationally very very hard i remember a project from years ago maybe 2012 or so we were working on a simulation for robotic arm placement in a constrained space we needed to find a valid configuration for the arm so its joints `theta1, theta2, theta3...` had to satisfy a bunch of inequalities some were simple like joint limits `-pi < theta1 < pi` but others were gnarly non-linear ones to avoid self-collision and obstacles involving sines and cosines all over the place `sqrt((x1-x2)^2 + (y1-y2)^2) > d_min` where `x1, y1` depended on `sin(theta1)` and `cos(theta2)` etc you get the picture it was a mess my first instinct back then as a younger dev was to build something myself i thought hey i'll just do some kind of random sampling or a clever gradient descent to push the arm away from invalid states i wasted three weeks literally three weeks building a custom solver that would get stuck in local minima all the time or just oscillate forever without finding a solution it was a nightmare my senior dev finally walked over saw what i was doing and basically laughed he told me to never ever roll your own optimizer unless you're a PhD in numerical analysis and even then you should think twice he pointed me to existing libraries and within a day i had a working prototype that was a hard lesson learned you don't build your own database for a production app you don't write your own crypto library and you definitely don't write your own non-linear constraint solver unless it's the entire point of your job so let's get to the practical stuff how you actually solve this today there are a few common approaches depending on your exact problem and what ecosystem you're in i'll assume you're probably using python since that's where a lot of this work happens these days the first and most common entry point for people is `scipy.optimize.minimize` it's a general purpose tool it's not strictly an inequality solver it's an *optimizer* but you can trick it into being a solver what you do is you define a dummy objective function like `f(x) = 0` and then you provide your inequalities as a set of constraints the solver's job is then to find a point that satisfies all the constraints while minimizing the dummy function which is trivial this effectively just becomes a search for a feasible point here's how you might do it let's say we want to find a point `(x, y)` that satisfies two conditions 1 `x^2 + y^2 < 25` (inside a circle of radius 5) 2 `y > sin(x)` (above the sine wave) you have to rewrite your inequalities into the form `g(x) >= 0` so they become 1 `25 - x^2 - y^2 > 0` 2 `y - sin(x) > 0` here's a code snippet using scipy it should just work if you have scipy installed ```python import numpy as np from scipy.optimize import minimize # We want to find a point (x, y) that satisfies: # 1. x^2 + y^2 < 25 --> 25 - (x^2 + y^2) > 0 # 2. y > sin(x) --> y - sin(x) > 0 # The constraints need to be in a list of dictionaries. # For 'ineq' type, the function must be non-negative (>= 0). cons = ( {'type': 'ineq', 'fun': lambda p: 25 - (p[0]**2 + p[1]**2)}, {'type': 'ineq', 'fun': lambda p: p[1] - np.sin(p[0])} ) # We don't care about optimizing anything, we just want a feasible point. # So we can just minimize a constant function. objective_func = lambda p: 0 # We need an initial guess to start the search. # Sometimes the choice of x0 matters a lot if the feasible region is weird. x0 = np.array([1.0, 1.0]) # Now we call the solver. SLSQP is a good algorithm for this. # Sequential Least Squares Programming result = minimize(objective_func, x0, method='SLSQP', constraints=cons) if result.success: found_point = result.x print(f"Successfully found a feasible point: x={found_point[0]:.4f}, y={found_point[1]:.4f}") # You can verify the constraints here val1 = 25 - (found_point[0]**2 + found_point[1]**2) val2 = found_point[1] - np.sin(found_point[0]) print(f"Constraint 1 value (must be > 0): {val1:.4f}") print(f"Constraint 2 value (must be > 0): {val2:.4f}") else: print("Solver failed to find a solution.") print(f"Message: {result.message}") ``` this is your bread and butter approach it's versatile but it might be slow and it can fail to find a solution even if one exists especially if the problem is highly non-convex or you provide a bad initial guess now if your problem is special if it's a *convex* problem then you can use much more powerful tools these tools are faster and they are guaranteed to find a solution if one exists convex optimization is a whole field but the short version is if your inequality `f(x) <= 0` has `f` as a convex function and your feasible set forms a convex shape then you're in business a library like `CVXPY` is fantastic for this it uses something called disciplined convex programming you basically state your problem in a natural way and cvxpy figures out if it's convex and if so it passes it to a powerful backend solver like SCS or ECOS let's reframe our problem slightly to be convex friendly maybe we want `x^2 + y^2 <= 25` and `exp(y) <= x` (this second one isn't but just for the sake of showing the syntax) CVXPY forces you to think about the problem structure Why did the programmer quit his job? Because he didn't get arrays. here's how you might write it in cvxpy notice how it's more declarative you're describing the problem not the steps to solve it ```python import cvxpy as cp import numpy as np # We want to find a point (x, y) that satisfies convex constraints # This is a bit of a contrived example to show syntax # 1. x^2 + y^2 <= 25 (This is a convex set, so quad_over_lin is a good way to express it) # 2. y >= 2*x + 1 (This is a linear and thus convex inequality) # Define the variables x = cp.Variable() y = cp.Variable() # Define the list of constraints constraints = [ cp.sum_squares(cp.hstack([x, y])) <= 25, # Expresses x^2+y^2 <= 25 y >= 2*x + 1 ] # Create a problem. We just want to find a feasible solution # so we use a dummy objective. prob = cp.Problem(cp.Minimize(0), constraints) # Solve the problem prob.solve() if prob.status == 'optimal': print(f"Successfully found a feasible point: x={x.value:.4f}, y={y.value:.4f}") # Verify print(f"x^2 + y^2 = {(x.value**2 + y.value**2):.4f} (should be <= 25)") print(f"y - (2x+1) = {(y.value - (2*x.value + 1)):.4f} (should be >= 0)") else: print(f"Could not solve the problem. Status: {prob.status}") ``` CVXPY is amazing when your problem fits its model if it doesn't fit the DCP ruleset it will complain and not let you run it which is actually a feature because it saves you from trying to solve an intractable problem with the wrong tool finally there's a third class of solvers for when your problem is less about continuous optimization and more about logic and satisfiability this is where you get into SMT solvers Satisfiability Modulo Theories the big name here is Z3 from Microsoft Research Z3 is a beast it can solve systems of constraints over real numbers integers bit-vectors all sorts of things it thinks about the problem not as finding a minimum of a function but as finding a model that makes a set of logical assertions true it's a different mindset let's take our original problem again and solve it with z3 ```python from z3 import Solver, Real, And, sat # We want to find a point (x, y) that satisfies: # 1. x^2 + y^2 < 25 # 2. y > sin(x) --> Z3 doesn't have sin directly, so this becomes an issue. # This shows a limitation of SMT solvers: they are for algebraic constraints. # Let's change the problem to something Z3 is good at: # 1. x^2 + y^2 < 25 # 2. y > x^2 - 4 # Create Z3 variables. These are symbolic variables. x = Real('x') y = Real('y') # Create a solver instance s = Solver() # Add constraints to the solver s.add(x**2 + y**2 < 25) s.add(y > x**2 - 4) # Check if the constraints are satisfiable result = s.check() if result == sat: # If it is satisfiable, get the model (a solution) m = s.model() # Z3 gives arbitrary-precision rational numbers, we can convert to float x_val = m[x].as_fraction() y_val = m[y].as_fraction() print(f"Successfully found a feasible point: x={float(x_val):.4f}, y={float(y_val):.4f}") else: print("The constraints are unsatisfiable (unsat) or the solver timed out (unknown).") ``` you can see Z3 is incredibly powerful for polynomial constraints and logic if you had `And(Or(x>5, y<2), z*x==1)` kind of stuff Z3 would be the tool to reach for its weakness is transcendental functions like `sin` `cos` `exp` which are the bread and butter of the numerical optimization world so my advice is this 1 start with `scipy.optimize.minimize` it's the general swiss army knife 2 if you can prove your problem is convex learn about `CVXPY` it will be much more reliable and faster 3 if your problem has a lot of logic integer variables or is more about satisfiability than optimization look into an SMT solver like `Z3` instead of just dropping links that will be dead in a year you should really get some foundational knowledge go find a copy of `Convex Optimization` by Boyd and Vandenberghe it's basically the bible for this stuff seriously every person in this field has read it for the hardcore numerical details the book series `Numerical Recipes` is a classic for a reason it explains the guts of algorithms like SLSQP and if you want to understand how these solvers actually work under the hood look up papers on `Interior-Point Methods` and `Sequential Quadratic Programming` or SQP that's the magic behind a lot of this don't try to write your own solver unless it's for a class project you'll just end up in a world of pain trust me on that one use the tools built by people who have dedicated their lives to this stuff good luck