How to generate constraints dynamically in scipy.optimize?

Question

Well what I was trying to do was to model the following using scipy.optimize.minimize.

What I'm trying to optimize is this function with its constraints: Here variable V is a list of variables, list's length is equal to the size of Omega. What I did so far is the following:

import numpy as np
from scipy.linalg import norm
from scipy.optimize import minimize

# set of concepts
M = ['linear algebra','seq2seq', 'artificial neural network','pointer networks']

#subchapters
S1=['linear algebra', 'differential calculus','backpropagation']
S2=['linear algebra','seq2seq', 'artificial neural network']
S3=['linear algebra','seq2seq', 'artificial neural network','pointer networks']

#vector representation of the subchapter in the concept space
x=[[1,1,0,0],
   [1,1,1,0],
   [1,1,1,1]]

# set of prerequisite relations among subchapters (entered manually for testing)
Omega = [(1, 2),(2,3),(1,3)]

# number of concepts
m = len(M)

# define of theta and lambda constants (manual)
theta = 2
lamda = 1

# matrix A is a m x m matrix , where m is the number of concepts
# it represents the prerequisite relations among concepts
# A is generated randomly
np.random.seed(43)
#A = np.zeros((m,m), dtype = int)
A = np.random.randint(2, size=(m,m))

# define the slack variable V as an array of size equal to len(Omega)
V = np.empty(len(Omega), dtype=float)

bnds=[]
# bounds -1 and 1 , create the array
# -1 <= a[i][j] <= 1
bnds_a_s_t = [bnds.append((-1, 1)) for _ in range(np.size(A))]
# bounds for the slack variable V, V is positive
bnds_V_i_j = [bnds.append((0,np.inf)) for _ in range(np.size(V))]

#constraints
cons=[]
#equality constraint
#a[s][t] + a[t][s] = 0
def equality_constraint(X):
    A_no_flatten=X[:len(X)-len(Omega)]
    #reconstruct matrix A
    A=np.reshape(A_no_flatten,(m,m))

    for s in range(m):
        for t in range(m):
            r=A[s][t]+A[t][s]
            #r=0 constraint
            con = {'type': 'eq', 'fun': lambda X: r} 
            cons.append(con)

# inequality constraint
#x[i].T @ (C[i][j] * A) @ x[j]
def inequality_constraint(X,x):
    for couple in Omega:
        # get the i and j
        i = couple[0]
        j = couple[1]
        
        #initialize C to 1s
        C = np.full((m,m), 1, dtype = int)
        # take all elements from X except last len(Omega) elements
        A_no_flatten=X[:len(X)-len(Omega)]
        # reconstruct list V
        V=X[-len(Omega):]
        #index for V
        f=0
        #reconstruct matrix A
        A=np.reshape(A_no_flatten,(m,m))
        
        #construct C[i][j]
        for s in range(m):
            for t in range(m):
                if x[i][t]>0 or x[j][s]>0:
                    C[s][t]=0
                else:
                    C[s][t]=1
        
        first= x[i].T
        second = C*A
        third = x[j]
        
        first_sec = first@second
        res=first_sec@third
        ineq_con = {'type': 'ineq', 'fun': lambda X: res -theta +V[f]}
        f+=1
        cons.append(ineq_con)
        
# arguments passed to the function
# here we pass x matrix 
# arguments are passed and used in constraints and in the objective function
# the objective function will minimize A and V which are matrix A and slack variable V
arguments=(x,)

# objective function
def objective(X, x):
    A_no_flatten=X[:len(X)-len(Omega)]
    # reconstruct list V
    V=X[-len(Omega):]
    #reconstruct matrix A
    A=np.reshape(A_no_flatten,(m,m))
    
    # sum of square V
    sum_square=0.0
    for it in V:
        sum_square+=it**2
    # sum of square V * lambda
    sum_square_lambda=sum_square*lamda
    
    return norm(A, 1) + sum_square_lambda

# list of variables to pass to the objective function
#pass the x0.flatten() which is the A + V combined, and then when in objective function we recreate them
# the first one A is all except the last s items where s is the size of V
# and then V is the rest

B = A.flatten()
p0 = np.append(B,V)

# scipy minimize
sol = minimize(objective, x0 = p0, args=arguments, bounds=bnds, constraints=cons)
print(sol.x)

What I get is the following:

[-7.73458566e-010  0.00000000e+000  4.99999997e-001  1.00000000e+000
  1.00000000e+000  0.00000000e+000 -5.00000003e-001  1.00000000e+000
  1.00000000e+000  1.00000000e+000  4.99999997e-001 -7.73458566e-010
 -7.73458566e-010  0.00000000e+000  4.99999997e-001 -7.73458566e-010
  6.01347002e-154  1.07176259e-311  0.00000000e+000]

Which doesn't respect the constraints and is not what I expected

What I don't know is that is it correct to add constraints like that, because I don't seem to call the constraints function, and I need to add them in a loop, and each function depends on X which is the list to minimize. When I print the cons array it is empty and I know, but I didn't find another way to add the constraint a[s][t]+a[t][s]=0 and the other one, I don't know if my approach is correct, thank you for your help in advance, much appreciated.

Answer 1

This isn't a complete answer, but it may get you started. As already mentioned in the comments, your list of constraints cons is empty when passed to the minimize method. So let's consider your first equality constraint. There are a few problems:

• Each time you call the function equality_constraint, you'd append new constraints to your list cons.

• Since you pass each constraint A[s][t] + A[t][s] == 0 as scalar function, it is quite cumbersome. Instead, you could use a single vector-valued function:

def constraint1(X):
    A_no_flatten = X[:len(X)-len(Omega)]
    A = np.reshape(A_no_flatten,(m,m))
    return A.flatten() + A.T.flatten()

As the name indicates, the .flatten() method flattens the matrix to a vector and A.T is just the transpose of A. Now you can add this constraint:

cons = []
cons.append({'type': 'eq', 'fun': constraint1})

Proceed similarly for the other constraint.