pythonmatrixoptimizationmathematical-optimizationgekko

Using GEKKO to optimize two matrix/vector equations


I wanted to use Gekko to solve an optimization (production mix) problem, I have a few numpy arrays, which I want to use in some vectorized equation.

They idea is, I have two simple matrix equations:

Ax < b

and

Cx = y

Where I will use previously prepared (const) numpy arrays for A,b,C.

Also, Y is a variable (int). X is a 1D array.

But if 'C' , 'X' , are 1D Arrays, and 'Y' is just a regular int variable (resulted from something like multiplying 1x100 matrix by 100x1 matrix), so I cant use m.axb for the second equation (It keeps giving me errors).

So I try to use them in a simple GEKKO script like this:

m = GEKKO(remote=False)

x = m.axb(const_machXparts, const_Mach_Max_Cap , etype='<')
y = m.Intermediate(const_FP_Wt_Matrix @ x)

m.Maximize(y)
m.solve(disp=True)
x.value

But I get the error: Exception: This steady-state IMODE only allows scalar values.

To my understanding, I don't need to define x and y using:

x = m.Array(m.Var, len(part_matrix_.columns), lb=0)

y = m.Var(lb=0)

Because I already have them defined using another way, I'm not sure if this is the cause of the problem, but adding these two lines to the beginning of the script, does not solve it.

I tried changing IMODE values using m.options.IMODE = 3,4,5 ..etc but it doesn't change anything.

I tried doing

x = m.Array(m.Var, len(part_matrix_.columns), lb=0)

then

m.axb(const_machXparts, const_Mach_Max_Cap , x=x , etype='<')

But it gives an "int have no length" error. Which only goes away when I assign the whole expression to x. Also documentation says: Usage: x = m.axb(A,b,etype='=,<,>,<=,>=',sparse=[True,False])

So now I'm stuck with the lines:

x = m.axb(const_machXparts, const_Mach_Max_Cap , etype='<')

y = m.Intermediate(const_FP_Wt_Matrix @ x)

I don't know what else to try.

Any help is appreciated. Thank you very much for your time in advance.


Solution

  • Here is an example with sample values for A, b, and C.

    import numpy as np
    from gekko import GEKKO
    
    # Initialize model
    m = GEKKO(remote=False)
    
    # Given matrices and constants
    A = np.array([[1, 2], [3, 4]])   # Example matrix A
    b = np.array([55, 41])           # Example vector b
    C = np.array([2, 3])             # Example row vector C
    
    # Constraint: Ax <= b
    x = m.axb(A,b,x=None,etype='<=',sparse=False)
    
    # alternative definition
    #Ar,Ac = A.shape
    #x = m.Array(m.Var, Ac, lb=0)  # Array of GEKKO variables for x (e.g., 1D array)
    #for i in range(A.shape[0]):
    #    m.Equation(m.sum([A[i, j] * x[j] for j in range(A.shape[1])]) <= b[i])
    
    # Equation: Cx = y
    yi = m.Var(lb=0,ub=10,integer=True)  # Scalar variable
    m.Equation(m.sum([C[i] * x[i] for i in range(len(C))]) == yi)
    
    # Objective: Maximize yi
    m.Maximize(yi)
    
    # Solve the problem
    m.options.SOLVER = 1
    m.solve(disp=True)
    
    # Results
    print('Optimized x:', [xi.value[0] for xi in x])
    print('Optimized y:', yi.value[0])
    

    I've also included an alternative definition (commented out) if you'd like to define the inequality constraint with a list comprehension. The second equation with C*x = y can't be used with the m.axb() function because y is a variable.

    The solution to this problem is:

     ----------------------------------------------------------------
     APMonitor, Version 1.0.3
     APMonitor Optimization Suite
     ----------------------------------------------------------------
    
    
     --------- APM Model Size ------------
     Each time step contains
       Objects      :  2
       Constants    :  0
       Variables    :  6
       Intermediates:  0
       Connections  :  5
       Equations    :  4
       Residuals    :  4
    
     Number of state variables:    6
     Number of total equations: -  6
     Number of slack variables: -  0
     ---------------------------------------
     Degrees of freedom       :    0
    
     ----------------------------------------------
     Steady State Optimization with APOPT Solver
     ----------------------------------------------
    Iter:     1 I:  0 Tm:      0.00 NLPi:    1 Dpth:    0 Lvs:    2 Obj: -3.86E-01 Gap:       NaN
    --Integer Solution:  -1.00E+01 Lowest Leaf:  -1.00E+01 Gap:   0.00E+00
    Iter:     2 I:  0 Tm:      0.00 NLPi:    2 Dpth:    1 Lvs:    2 Obj: -1.00E+01 Gap:  0.00E+00
     Successful solution
    
     ---------------------------------------------------
     Solver         :  APOPT (v1.0)
     Solution time  :  0.019 sec
     Objective      :  -10.
     Successful solution
     ---------------------------------------------------
    
    
    Optimized x: [2.3529411765, 1.7647058824]
    Optimized y: 10.0