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
visithttps://groups.google.com/d/msgid/sympy/403209cd-1bf6-4284-aef8-bc3a313a23e
bn%40googlegroups.com.
