3

I have a system of equations to solve which goes like this:

https://i.stack.imgur.com/FQY2g.jpg

the left side of the equation is the sum, and the right side is the constraints. I have a system of N nonlinear equations with N variables.

I tried solving it with scipy fsolve but it won't converge, sympy also has a solver, doesn't converge either. I came across pyscipopt which seem to work but behaves inconsistently and breaks down.


import numpy as np   
import networkx as nx
''' set up the problem with n nodes, and the constraint array c with the degree sequence'''

''' example : power law degree distribution''' 
n=8
c = np.round(nx.utils.random_sequence.powerlaw_sequence(n))


''' other examples '''
#n=8
#c=[1, 2, 2, 1, 2, 1, 1, 2]
#n=3
#c=[2,1,1]

from pyscipopt import Model, quicksum

m = Model()
''' tolerance for convergence'''
tol = 1e-6

print('create the lagrange multipliers')
X = dict(zip(range(n), [m.addVar(vtype="C", lb=0) for i in range(n)]))
#X = dict(zip(range(n), [m.addVar(vtype="C") for i in range(n)]))

'''  create the variable for the objective function because it's nonlinear
and must be linearized'''
Y = m.addVar(vtype='C')

print('set up the constraints')

''' the equations are essentially sum_i - degree_i <= tolerance'''

'''  set up the constraints as sums of the lagrange multipliers as written in the papers'''
system =[]
for i in range(n):
    idx = np.arange(n)
    idx = idx[idx!=i]
    k= quicksum(X[i]*X[j]/(1+X[i]*X[j]) for j in idx) -c[i]
    ''' the equations are essentially sum_i - degree_i <= tolerance'''
    m.addCons(k<=tol)
    system.append(k)

m.addCons( Y <= quicksum(system[j] for j in range(n)) )
m.setObjective(Y, 'minimize')
print('all constraints added, now optimization')
m.optimize()

So I have a system array where I store all the constraints to be optimized, the k values are the constraints, I perform a double loop over the entire data set.

1) are there better, or faster ways of achieving this? I'm guessing that for large N it wouldn't be efficient.

one example which works fast:

n=8 c=[1, 2, 2, 1, 2, 1, 1, 2]

but other examples, especially once I increase N (I use little n in the code), get stuck somehow in limbo.

edit: regarding the examples I have to say that any example should work. The hard constraint in designing an example is simply k_i <= N.

Lionel Yelibi
  • 43
  • 6
  • 1
    Could you please share the entire code? `X` and `Y` are not initialized. For what sizes of `n` do you experience issues? – mattmilten Mar 31 '19 at 12:52
  • 1. I edited the post and added the entire code 2. I think even n=8 with slightly different constraints (power law, see code) isn't producing results. I know the simple example ( n=8 c=[1, 2, 2, 1, 2, 1, 1, 2] ) works but I don't know why the others don't – Lionel Yelibi Mar 31 '19 at 17:28
  • 1
    This seems to be a tough problem. Especially, it's a purely continuous non-linear problem, so SCIP is probably not the best suited solver. Did you check something like BARON, ANTIGONE, Couenne or LINDO? – mattmilten Apr 01 '19 at 18:13
  • ha! I will try these then. Thanks. I will report back in a week or so ^^ – Lionel Yelibi Apr 02 '19 at 01:55
  • 1
    In general, the absence of upper bounds can make the problem even harder and can lead to very long solving times. – mattmilten Apr 02 '19 at 12:34

0 Answers0