Блог

У меня есть система уравнений для решения, которая выглядит следующим образом:

https://imgur.com/oxVdM10

левая часть уравнения — это сумма, а правая часть — ограничения. У меня есть система из N нелинейных уравнений с N переменными.

Я пытался решить это с помощью scipy fsolve, но это не сходится, у sympy также есть решатель, тоже не сходится. Я наткнулся на pyscipopt, который, похоже, работает, но ведет себя непоследовательно и выходит из строя.

 
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()

  

Итак, у меня есть системный массив, где я храню все ограничения, подлежащие оптимизации, значения k являются ограничениями, я выполняю двойной цикл по всему набору данных.

1) существуют ли лучшие или более быстрые способы достижения этого? Я предполагаю, что для большого N это было бы неэффективно.

один из примеров, который работает быстро:

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

но другие примеры, особенно когда я увеличиваю N (я использую мало n в коде), каким-то образом застревают в неопределенности.

редактировать: что касается примеров, я должен сказать, что любой пример должен работать. Жестким ограничением при разработке примера является просто k_i <= N.