0

I have a nonlinear optimization problem which makes use of 3 decision variables, one of these variables is a single number (t), one is a vector with index i (S_i) and one is a matrix (Q_i,j) with indices i and j. I'm currently trying to use scipy.optimize.minimize to model and solve my problem but I can't get it to work.

Because of the multiple indices, constraints must often hold for all values of some index. My actual model is very large so let's consider this example which we assume is nonlinear:

Decision variable: Q_i,j

Objective:

minimize Sum over all i and all j, Q_i,j

Constraint_1:

Q_i,j / 2 >= 10 for all i, for all j

Current code I try to use:

from scipy.optimize import minimize
import numpy as np

I = 5
J = 5
x0 = np.zeros(I*J)

def obj(Q_ij):
   Q_ijcp = np.reshape(Q_ijcp,(I,J))
   return sum(Q_ij[i,j] for i in range(I) for j in range(J))

def cons_1(Q_ij):
   Q_ijcp = np.reshape(Q_ijcp,(I,J))
   return (Q_ij[i,j] / 2 - 10 for i in range(I) for j in range(J))

b = (0, 100)
bounds = []
for i in range(I):
   for j in range(J):
      bounds.append(b)

constraint1 = {"type": "ineq", "fun": cons_1}
constraints = [constraint1]

solution = minimize(obj, x0=x0, bounds=bounds, constraints=constraints, method='SLSQP')

Based on the user guide I found that for each constraint one must make a definition such that it can be entered into the solver which I try above but it does not work, how can I model this such that I don't have to make a definition for each possible i and j value in Constraint_1? (Such that I don't end up with 5 * 5 constraints as I=J=5)

Or are there any other packages with good documentation/examples which are easier to use in the case of using vectors and matrices with constraints for their indices?

Bas R
  • 175
  • 7
  • There's nothing in your problem implying that Q needs to be 2x2. You can treat it as a simple raveled array. `obj` is just `Q.sum`, and `cons_1` is something like `(Q > 20).sum()` or `((Q - 10) / 2).sum()` – Mad Physicist Dec 24 '21 at 20:39
  • In `cons_1` I do not wish to sum over the matrix _Q_, the constraint must hold for each value in the matrix _Q_ which is indexed by _i_ and _j_. So for an arbitrary _i_ and _j_ the corresponding value of _Q[i,j]_ divided by 2 must be larger than or equal to 10. For instance, _i = 2_ and _j = 3_, we assume _Q[2,3] = 30_, _30 / 2 = 15_ and thus the constraint would hold for _Q[2,3]_. – Bas R Dec 24 '21 at 21:18
  • How is contraint actually documented to work? Does it return False if any of the values is out of bounds? – Mad Physicist Dec 24 '21 at 21:22
  • No, perhaps my example wasn't very clear. The values for the matrix _Q_ must be chosen such that 1. the sum of them all is minimized and 2. all of the values idividually divided by 2 must be at least 10, i.e. each value must be at least 20. So this optimization problem would set all values in the matrix _Q_ to 10. My actual problem is ofcourse different than this, but it follows a similar pattern. – Bas R Dec 24 '21 at 21:56
  • I completely understood that part. What I'm asking you is to give some indication of what the documentation says `const_1` should return. Is it a float? a boolean? etc. – Mad Physicist Dec 24 '21 at 21:58
  • Sorry for the misunderstanding, I believe the expected returned value is just a float. – Bas R Dec 24 '21 at 22:38
  • You're almost there. A float representing what exactly? What are the expected values? Please read the docs carefully. – Mad Physicist Dec 24 '21 at 22:39
  • Sorry, I'm quite new to this and I don't fully get the question. But I believe the expected return is the outcome of the constraint equation based on some value of Q[i,j] which the model uses to ensure the constraint gets satisfied when it's searching for the optimal values to fill matrix Q. – Bas R Dec 24 '21 at 22:50
  • https://docs.scipy.org/doc/scipy/tutorial/optimize.html#sequential-least-squares-programming-slsqp-algorithm-method-slsqp these are the docs, perhaps it helps in answering your question as I do not fully understand it. – Bas R Dec 24 '21 at 22:52
  • "My actual model is very large" Scipy's nonlinear solvers are not really designed for large, sparse problems. – Erwin Kalvelagen Dec 26 '21 at 08:00

1 Answers1

0

Or are there any other packages with good documentation/examples which are easier to use in the case of using vectors and matrices with constraints for their indices?

Below is an example problem with Python gekko with solver IPOPT to solve a similar nonlinear optimization problem.

from gekko import GEKKO
m = GEKKO()
n = 5
Q = m.Array(m.Var,(n,n),lb=10,ub=100)
t = m.Var(lb=0,ub=1)
S = m.Array(m.Var,n,lb=0,ub=1)

# add constraint t*Q*S > 100
QS = Q@S # dot product
m.Equations([t*QS[i]>100 for i in range(n)])

m.Minimize(m.sum([Q[i,j] for i in range(n)
                         for j in range(n)]))
m.options.SOLVER=3
m.solve()

print(t.value[0])
print(S)
print(Q)

Optimizers APOPT and IPOPT can handle 100,000+ variables. Gekko Documentation and examples are available.

John Hedengren
  • 12,068
  • 1
  • 21
  • 25