That system is helpful and you've got some good advice already. I would add
that this sort of system can be reduced to a linear programming form and
sympy has a function that handles your system:
```
from sympy import *
from sympy.solvers.simplex import lpmin
def lpsys(eqs):
eqs = list(eqs)
maxes = Tuple(*eqs).atoms(Max)
syms_max = symbols(f'm0:{len(maxes)}', real=True)
rep_m = dict(zip(maxes, syms_max))
eqs2 = [e.xreplace(rep_m) for e in eqs] # linear now after replacing
Max with symbols
max_eq=[Eq(*z) for z in zip(syms_max, maxes)]
return eqs2, max_eq
eqs2, max_eq =lpsys(eqs) # your list of 239 equations is eqs
obj = sum(ms.lhs for ms in max_eq) # minimize the sum of m variables used
to represent the Max functions
cons = list(eqs2)
# m_i >= each argument for Max's
for meq in max_eq:
for a in meq.rhs.args:
cons.append(meq.lhs >= a)
opt_val, sol = lpmin(obj, cons)
```
This gives 22.9893360923961 for the opt_val and solutions for the 269
equations (30 maxes and 239 others).
/c
On Sunday, August 24, 2025 at 9:30:38 AM UTC-5 Oscar wrote:
> On Sun, 24 Aug 2025 at 12:01, Gábor Horváth <[email protected]> wrote:
> >
> >
> > > Could you post an example of a large system?
> >
> > Attached one with 239 equations to this email. Not sure if attached files
> > go through, though. If not, I can paste it directly into the email body,
> > it's size is 17K.
>
> Yes, it has come through.
>
> Okay looking at this you definitely want a special algorithm for these
> kinds of systems rather than just using a general purpose solver.
> Ideally SymPy would have a reduce function and if it did we could
> consider whether it should be able to handle larger systems of this
> kind efficiently.
>
> First find the Max expressions and any symbols that are in them:
>
> maxes = Tuple(*eqs).atoms(Max)
> special = Tuple(*maxes).free_symbols
>
> Then replace all Max expressions with symbols themselves:
>
> syms_max = symbols(f'm:{len(maxes)}')
> rep_m = dict(zip(maxes, syms_max))
> eqs2 = [eq.xreplace(rep_m) for eq in eqs]
>
> Now you have a linear system of equations for the original variables
> v0, v1 etc and also new variables m0, m1 representing the maxes. We
> want to get the smallest possible system involving only the maxes, the
> symbols that are in the maxes (special) and as few of the other
> symbols as possible.
>
> Choose the ordering of the symbols carefully and ask linsolve to
> reduce the underdetermined system:
>
> dontcare = [s for s in syms if s not in special]
> unknowns = list(syms_max) + dontcare + list(special)
> [sol] = linsolve(eqs2, unknowns)
>
> In so far as possible this tries to solve for syms_max in terms of
> everything else and to solve for everything in terms of the symbols in
> special. It succeeds in expressing all of the maxes in terms of
> everything else:
>
> >>> sol.atoms(Symbol) - set(syms)
> set()
>
> Mostly this expresses everything in terms of the special symbols
> although a few others still appear as free parameters:
>
> >>> print(Tuple(*sol).free_symbols)
> {v232, v238, v121, v223, v183, v197, v62, v3, v218, v189, v106, v214,
> v216, v231, v11, v213, v6, v88, v71, v120, v111, v180, v21, v129,
> v208, v44, v153, v192, v201, v184}
>
> >>> print(Tuple(*sol).free_symbols - special)
> {v238, v197, v129, v223, v192, v214}
>
> Everything is now expressed in terms of these 30 symbols of which 24
> are symbols appearing in the maxes and 6 are not. These 6 symbols
> appearing in maxes have been solved purely from the linear equations
> but are given in terms of the 30 symbols left over:
>
> >>> print(special - Tuple(*sol).free_symbols)
> {v25, v9, v154, v160, v229, v84}
>
> I think at this point then you can translate the solution into 30
> equations for the remaining 30 unknowns:
>
> rep = dict(zip(unknowns[30:], sol[30:]))
> eqs3 = [Eq(m.xreplace(rep), s) for s, m in zip(sol[:30], maxes)]
>
> Now you have 30 equations for 30 unknowns involving 30 maxes and the
> other equations and variables are all eliminated and solved
>
> The other equations after these 30 can be inspected to give bounds on
> the 30 variables left over. You have a parametric solution for the 209
> other variables in terms of the 30 symbols. You can substitute that
> solution into any auxiliary inequality constraints to further
> constraint the values of the 30 symbols (or prove unsatisfiability).
>
> These are the remaining Maxes:
>
> In [222]: for m in Tuple(*eqs3).atoms(Max): print(m.n(3))
> Max(0.802, v111 - v153 + v180)
> Max(0.757, v62)
> Max(0.597, 49.8*v153 - 30.8*v180 - 54.3*v183 + 210.0*v189 - 210.0*v201
> + 15.4*v21 - 130.0*v218 - 4.23*v3 - 21.7*v44 + 210.0*v88 - 33.9)
> Max(0.879, v189 - v201 + v88)
> Max(0.732, v21)
> Max(0.735, v216)
> Max(0.767, v44)
> Max(0.645, v231)
> Max(0.721, v184)
> Max(0.789, v121)
> Max(0.897, v189)
> Max(0.741, v3)
> Max(0.835, v153)
> Max(0.827, v88)
> Max(0.677, -100.0*v153 + 164.0*v180 + 46.8*v183 - 183.0*v189 +
> 183.0*v201 - 46.9*v21 + 128.0*v218 + 2.33*v3 - 11.2*v44 - 183.0*v88 +
> 1.07)
> Max(0.846, v201)
> Max(0.786, v183)
> Max(0.783, v6)
> Max(0.727, v213)
> Max(0.719, v71)
> Max(0.78, v111 - v153 + v218)
> Max(0.658, v208)
> Max(0.804, 2.92*v153 - 6.78*v180 - 0.734*v183 + 6.25*v189 - 1.72*v201
> + 1.16*v21 - 1.2*v218 - 0.0219*v3 - 0.123*v44 + 1.72*v88 - 0.473)
> Max(0.703, v11)
> Max(0.786, v120)
> Max(0.683, v232)
> Max(0.839, v218)
> Max(0.837, v106)
> Max(0.862, v180)
> Max(0.777, v111)
>
> These are given linearly in terms of the same symbols e.g.:
>
> >>> print(eqs3[0].n(3))
> Eq(Max(0.757, v62), 474.0*v106 + 0.503*v111 - 81.4*v120 - 312.0*v121
> - 953.0*v153 + 1.8e+3*v180 + 506.0*v183 - 2.1e+3*v189 + 2.03e+3*v201 -
> 1.91*v208 - 543.0*v21 - 0.416*v216 + 1.53e+3*v218 + 7.85*v232 -
> 1.58*v238 - 71.3*v3 - 213.0*v44 - 13.0*v6 + 28.1*v62 - 2.03e+3*v88 -
> 61.2)
>
> Maybe we are actually putting this the wrong way round. Keep the m
> symbols and invert to solve for the remaining symbols in terms of the
> Maxes:
>
> eqs4 = [Eq(m, s) for s, m in zip(sol[:30], syms_max)]
> [sol2] = linsolve(eqs4, syms2)
>
> Yes, that's better. Now you have solved for 209 symbols in terms of
> these 30 and you have solutions for these 30 in terms of the maxes.
> Now each time you pick a branch for a Max it directly tells you the
> value of some symbol which you can then eliminate from everything
> else. There is no further need to call linsolve because the linear
> equations are all solved and it is just a case of considering the
> maxes in turn.
>
> Or something like that...
>
> I guess this is already where you are at and you just need a more
> efficient way of pruning the 2^30 branches.
>
> --
> Oscar
>
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion visit
https://groups.google.com/d/msgid/sympy/403209cd-1bf6-4284-aef8-bc3a313a23ebn%40googlegroups.com.
