Gekko: Resource optimisation with allocation matrix and weights array

Question

I am working on a format of resource optimization problem and writing a GEKKO code to solve. Problem statement is as follows: Let assume there are 2 workers and 4 tasks. Each worker gets some rewards for the tasks it picks. Each task can be assigned to only one worker and each worker can only pick one tasks. The objective is to maximize the total rewards obtained by both workers.

As per the requirement of problem decision variables are as follows:

1. An 0-1 allocation matrix that decides which task to be assigned to which worker
2. A weights array, that assign relative weights to each task.

Reward given to worker is calculated based on the ratio of weights of task picked and sum of weights of all tasks allocated to worker.

Following is the code for same:

import numpy as np
from gekko import GEKKO

rewards = np.array([[2,5,8,10],[1,5,7,11]])

m = GEKKO(remote=False)

allocation = m.Array(m.Var,(2,4),lb=0,ub=1, integer=True)
weights = m.Array(m.Var,4,lb=0,ub=1)


def reward(allocation, weights, rewards):
    temp=np.copy(allocation)
    #sum_ = np.sum(weights)
    for i in range(allocation.shape[1]):
        temp[:, i] *= weights[i]#/sum_
    
    total_rewards = np.sum(temp.flatten() * rewards.flatten())
    return total_rewards


for j in range(4):
    m.Equation(m.sum(allocation[:,j])<=1) 


m.Maximize(reward(allocation, weights, rewards))

m.options.SOLVER = 1        # change solver (1=APOPT,3=IPOPT)
#m.open_folder()
m.solve()

print('allocation', allocation)
print('weights', weights)
print('Objective: ' + str(m.options.objfcnval))

In the reward function, I have commented "sum_" variable. In this case the output is:

----------------------------------------------------------------
 APMonitor, Version 1.0.1
 APMonitor Optimization Suite

 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            0
   Variables    :           16
   Intermediates:            0

   Connections  :            0
   Equations    :            5
   Residuals    :            5
 
 Number of state variables:             16
 Number of total equations: -            4
 Number of slack variables: -            4
 ---------------------------------------
 Degrees of freedom       :              8
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.00 NLPi:    4 Dpth:    0 Lvs:    3 Obj: -1.63E+01 Gap:       NaN
--Integer Solution:   0.00E+00 Lowest Leaf:  -1.63E+01 Gap:   2.00E+00
Iter:     2 I:  0 Tm:      0.00 NLPi:    2 Dpth:    1 Lvs:    2 Obj:  0.00E+00 Gap:  

2.00E+00
Iter:     3 I:  0 Tm:      0.00 NLPi:    3 Dpth:    1 Lvs:    4 Obj: -2.60E+01 Gap:  2.00E+00
--Integer Solution:  -2.60E+01 Lowest Leaf:  -2.60E+01 Gap:   0.00E+00
Iter:     4 I:  0 Tm:      0.00 NLPi:    1 Dpth:    2 Lvs:    4 Obj: -2.60E+01 Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)

 Solution time  :   1.729999999224674E-002 sec
 Objective      :   -26.0000000000000     
 Successful solution
 ---------------------------------------------------
 

allocation [[[1.0] [1.0] [1.0] [0.0]]
 [[0.0] [0.0] [0.0] [1.0]]]
weights [[1.0] [1.0] [1.0] [1.0]]
Objective: -26.0

But if i uncomment "sum_" variable i.e use ratio of weights and total weights of task allocated to worker. I am getting objective value as nan

----------------------------------------------------------------


APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            0
   Variables    :           16
   Intermediates:            0
   Connections  :            0
   Equations    :            5
   Residuals    :            5
 
 Number of state variables:             16
 Number of total equations: -            4
 Number of slack variables: -            4
 ---------------------------------------
 Degrees of freedom       :              8
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.01 NLPi:    2 Dpth:    0 Lvs:    0 Obj:       NaN Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)

     Solution time  :   2.390000000013970E-002 sec
     Objective      :                      NaN
     Successful solution
     ---------------------------------------------------
     
    
    allocation [[[0.0] [0.0] [0.0] [0.0]]
     [[0.0] [0.0] [0.0] [0.0]]]
    weights [[0.0] [0.0] [0.0] [0.0]]
    Objective: nan

Please let me know how to fix this issue or any pointer to debug this.

Thanks

Answer 1

Set a lower bound on sum_ to avoid divide-by-zero in the denominator.

sum_ = m.Var(lb=1e-3)
m.Equation(sum_ == np.sum(weights))

The NaN is because the default guess value is 0 and the summation is initially equal to zero as well. Here is a complete script that solves successfully:

import numpy as np
from gekko import GEKKO

rewards = np.array([[2,5,8,10],[1,5,7,11]])
m = GEKKO(remote=False)
allocation = m.Array(m.Var,(2,4),lb=0,ub=1, integer=True)
weights = m.Array(m.Var,4,lb=0,ub=1)

def reward(allocation, weights, rewards):
    temp=np.copy(allocation)
    sum_ = m.Var(lb=1e-3)
    m.Equation(sum_ == np.sum(weights))
    for i in range(allocation.shape[1]):
        temp[:, i] *= weights[i]/sum_    
    total_rewards = np.sum(temp.flatten() * rewards.flatten())
    return total_rewards

for j in range(4):
    m.Equation(m.sum(allocation[:,j])<=1) 
m.Maximize(reward(allocation, weights, rewards))

m.options.SOLVER = 1        # change solver (1=APOPT,3=IPOPT)
m.solve()

print('allocation', allocation)
print('weights', weights)
print('Objective: ' + str(m.options.objfcnval))

Here is the output:

 ----------------------------------------------------------------
 APMonitor, Version 1.0.0
 APMonitor Optimization Suite
 ----------------------------------------------------------------

 --------- APM Model Size ------------
 Each time step contains
   Objects      :  0
   Constants    :  0
   Variables    :  17
   Intermediates:  0
   Connections  :  0
   Equations    :  6
   Residuals    :  6

 Number of state variables:    17
 Number of total equations: -  5
 Number of slack variables: -  4
 ---------------------------------------
 Degrees of freedom       :    8

 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.00 NLPi:    7 Dpth:    0 Lvs:    3 Obj: -1.10E+01 Gap:       NaN
--Integer Solution:  -1.10E+01 Lowest Leaf:  -1.10E+01 Gap:   0.00E+00
Iter:     2 I:  0 Tm:      0.00 NLPi:    2 Dpth:    1 Lvs:    3 Obj: -1.10E+01 Gap:  0.00E+00
 Successful solution

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.02 sec
 Objective      :  -11.
 Successful solution
 ---------------------------------------------------

allocation [[[0.0] [0.0] [1.0] [0.0]]
 [[0.0] [0.0] [0.0] [1.0]]]
weights [[0.0] [0.0] [0.0] [0.0020037011152]]
Objective: -11.0

Answer 2

Based on the following:

• in the allocation matrix, task columns sum to 0 or 1;
• that means each element of the weight vector is assigned to 0 or 1 workers;
• the sum of assigned weights per worker row is always 1:

that means that the sum of all weights must always equal the number of workers, 2 in this case. You can successfully linearise with big-M constraints and auxiliary variables for weights as assigned to each worker-task pair. This looks like:

import numpy as np
import pulp

rewards = np.array((
    (2, 5, 8, 10),
    (1, 5, 7, 11),
))
workers = range(2)
tasks = range(4)

allocation = pulp.LpVariable.matrix(
    name='alloc', cat=pulp.LpBinary,
    indices=(workers, tasks),
)
weights = pulp.LpVariable.matrix(
    name='weight', cat=pulp.LpContinuous, lowBound=0, upBound=1,
    indices=tasks,
)
assigned_weights = pulp.LpVariable.matrix(
    name='aweight', cat=pulp.LpContinuous, lowBound=0,
    indices=(workers, tasks),
)

'''
In the allocation matrix, columns sum to 0 or 1.
That means each element of the weight vector is assigned to 0 or 1 workers.
That in turn means that the sum of all weights is equal to the number of workers.
'''
weight_sum = len(workers)

m = pulp.LpProblem(name='workload_distribution', sense=pulp.LpMaximize)
m.setObjective(
    pulp.lpSum(
        pulp.lpDot(weight_row, reward_row)/weight_sum
        for weight_row, reward_row in zip(assigned_weights, rewards)
    )
)

# Each task can be assigned up to once
for task in tasks:
    m.addConstraint(
        name=f'excl_t{task}',
        constraint=pulp.lpSum(
            row[task] for row in allocation
        ) <= 1,
    )

M = 10
for worker, alloc_row, weight_row in zip(workers, allocation, assigned_weights):
    # Per worker, assigned weights need to sum to 1
    m.addConstraint(
        name=f'weights_w{worker}',
        constraint=pulp.lpSum(weight_row) == 1,
    )

    # Weight assignment needs to track task assignment
    for task, alloc, aweight, weight in zip(tasks, alloc_row, weight_row, weights):
        suffix = f'_w{worker}_t{task}'
        m.addConstraint(
            name='weight_alloc' + suffix,
            constraint=aweight <= alloc,
        )
        m.addConstraint(
            name=f'weightweight_lo' + suffix,
            constraint=aweight <= weight,
        )
        m.addConstraint(
            name=f'weightweight_hi' + suffix,
            constraint=aweight >= weight - M*(1 - alloc),
        )

print(m)
m.solve()
assert m.status == pulp.LpStatusOptimal

print('Allocations:')
for row in allocation:
    print([int(a.value()) for a in row])

print('\nWeights:')
print([w.value() for w in weights])

print('\nAssigned weights:')
for row in assigned_weights:
    print([w.value() for w in row])

workload_distribution:
MAXIMIZE
1.0*aweight_0_0 + 2.5*aweight_0_1 + 4.0*aweight_0_2 + 5.0*aweight_0_3 + 0.5*aweight_1_0 + 2.5*aweight_1_1 + 3.5*aweight_1_2 + 5.5*aweight_1_3 + 0.0
SUBJECT TO
excl_t0: alloc_0_0 + alloc_1_0 <= 1

excl_t1: alloc_0_1 + alloc_1_1 <= 1

excl_t2: alloc_0_2 + alloc_1_2 <= 1

excl_t3: alloc_0_3 + alloc_1_3 <= 1

weights_w0: aweight_0_0 + aweight_0_1 + aweight_0_2 + aweight_0_3 = 1

weight_alloc_w0_t0: - alloc_0_0 + aweight_0_0 <= 0

weightweight_lo_w0_t0: aweight_0_0 - weight_0 <= 0

weightweight_hi_w0_t0: - 10 alloc_0_0 + aweight_0_0 - weight_0 >= -10

weight_alloc_w0_t1: - alloc_0_1 + aweight_0_1 <= 0

weightweight_lo_w0_t1: aweight_0_1 - weight_1 <= 0

weightweight_hi_w0_t1: - 10 alloc_0_1 + aweight_0_1 - weight_1 >= -10

weight_alloc_w0_t2: - alloc_0_2 + aweight_0_2 <= 0

weightweight_lo_w0_t2: aweight_0_2 - weight_2 <= 0

weightweight_hi_w0_t2: - 10 alloc_0_2 + aweight_0_2 - weight_2 >= -10

weight_alloc_w0_t3: - alloc_0_3 + aweight_0_3 <= 0

weightweight_lo_w0_t3: aweight_0_3 - weight_3 <= 0

weightweight_hi_w0_t3: - 10 alloc_0_3 + aweight_0_3 - weight_3 >= -10

weights_w1: aweight_1_0 + aweight_1_1 + aweight_1_2 + aweight_1_3 = 1

weight_alloc_w1_t0: - alloc_1_0 + aweight_1_0 <= 0

weightweight_lo_w1_t0: aweight_1_0 - weight_0 <= 0

weightweight_hi_w1_t0: - 10 alloc_1_0 + aweight_1_0 - weight_0 >= -10

weight_alloc_w1_t1: - alloc_1_1 + aweight_1_1 <= 0

weightweight_lo_w1_t1: aweight_1_1 - weight_1 <= 0

weightweight_hi_w1_t1: - 10 alloc_1_1 + aweight_1_1 - weight_1 >= -10

weight_alloc_w1_t2: - alloc_1_2 + aweight_1_2 <= 0

weightweight_lo_w1_t2: aweight_1_2 - weight_2 <= 0

weightweight_hi_w1_t2: - 10 alloc_1_2 + aweight_1_2 - weight_2 >= -10

weight_alloc_w1_t3: - alloc_1_3 + aweight_1_3 <= 0

weightweight_lo_w1_t3: aweight_1_3 - weight_3 <= 0

weightweight_hi_w1_t3: - 10 alloc_1_3 + aweight_1_3 - weight_3 >= -10

VARIABLES
0 <= alloc_0_0 <= 1 Integer
0 <= alloc_0_1 <= 1 Integer
0 <= alloc_0_2 <= 1 Integer
0 <= alloc_0_3 <= 1 Integer
0 <= alloc_1_0 <= 1 Integer
0 <= alloc_1_1 <= 1 Integer
0 <= alloc_1_2 <= 1 Integer
0 <= alloc_1_3 <= 1 Integer
aweight_0_0 Continuous
aweight_0_1 Continuous
aweight_0_2 Continuous
aweight_0_3 Continuous
aweight_1_0 Continuous
aweight_1_1 Continuous
aweight_1_2 Continuous
aweight_1_3 Continuous
weight_0 <= 1 Continuous
weight_1 <= 1 Continuous
weight_2 <= 1 Continuous
weight_3 <= 1 Continuous

...


Result - Optimal solution found

Objective value:                9.50000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01

Allocations:
[0, 0, 1, 0]
[0, 0, 0, 1]

Weights:
[0.0, 0.0, 1.0, 1.0]

Assigned weights:
[0.0, 0.0, 1.0, 0.0]
[0.0, 0.0, 0.0, 1.0]

It makes perfect sense that the best solution moves all of the weight to the highest reward elements of 8 and 11 respectively, with no weight allocated to the lower-reward elements. In fact, since you've said that "each worker can only pick one task", this further reduces to a weight-less assignment problem:

import numpy as np
import pulp

rewards = np.array((
    (2, 5, 8, 10),
    (1, 5, 7, 11),
))
workers = range(2)
tasks = range(4)

allocation = pulp.LpVariable.matrix(
    name='alloc', cat=pulp.LpBinary,
    indices=(workers, tasks),
)

m = pulp.LpProblem(name='workload_distribution', sense=pulp.LpMaximize)
m.setObjective(
    pulp.lpSum(
        pulp.lpDot(alloc_row, reward_row)
        for alloc_row, reward_row in zip(allocation, rewards)
    )
)

# Each task can be assigned to only one worker
for task in tasks:
    m.addConstraint(
        name=f'excl_t{task}',
        constraint=pulp.lpSum(
            row[task] for row in allocation
        ) <= 1,
    )

# Each worker can only pick one task
for worker, alloc_row in zip(workers, allocation):
    m.addConstraint(
        name=f'excl_w{worker}',
        constraint=pulp.lpSum(alloc_row) == 1,
    )

print(m)
m.solve()
assert m.status == pulp.LpStatusOptimal

print('Allocations:')
for row in allocation:
    print([int(a.value()) for a in row])

Result - Optimal solution found

Objective value:                19.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

Allocations:
[0, 0, 1, 0]
[0, 0, 0, 1]