2

I am working on a optimisation problem and I am using Gekko to solve it. Consider the scenario where there are 2 machines and 5 users and each machines has 2 resources to distribute to users. The user gain some points based on the number of resource allocated. Also there are certain demands from the user before resource is allocated. My objective is two fold, I want to maximise a combination of total gain and fairness of the system.

Fairness of user is defined as total gain divided by demand. I have two decision variable, one is allocation matrix between machines and user & other is weights vector for users. The resources are allocated in proportional to the weights of user. (The two decision variables are compulsory in the problem I am working on)

Following is the code i have worked upon:

import numpy as np
from gekko import GEKKO

def allocation(rewards, num_machines, num_users, num_resources, C1, C2, demands, cumulative_gain):
    
    
    m = GEKKO(remote=False)

    # Decision variables
    link = m.Array(m.Var, (num_machines, num_users), lb=0, ub=1, integer=True)
    weights = m.Array(m.Var, num_users, lb=0.1, ub=1)

    # Constraints
    for u in range(num_users):
        m.Equation(m.sum(link[:, u]) <= 1)  # each user is assigned to only 1 machine

    for b in range(num_machines):
        m.Equation(m.sum(link[b, :] * weights) <= 1)  # sum of weights of users allocated to each machine should be 1

    # Compute Fairness
    fairness, reward = compute_fairness(link, weights, rewards, demands, cumulative_gain, num_resources, m)

    # Objective
    m.Maximize(C1 * reward + C2 * fairness)

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

    optimized_weights = np.array([weight.value[0] for weight in weights])
    optimized_link = np.ndarray((num_machines, num_users))
    for i in range(link.shape[0]):
        for j in range(link.shape[1]):
            optimized_link[i][j] = link[i][j][0]

    return optimized_link, optimized_weights, m

def compute_fairness(link, weights, rewards, demands, cumulative_gain, num_resources, m):
    """
    Computes the fairness of the system based on the resources allocated to users using GEKKO-compatible logic.
    """
    num_machines = link.shape[0]
    num_users = link.shape[1]
    
    total_gain_u = []
    user_fairness = []
    total_machine_weight = []
    
    for b in range(num_machines):
        total_machine_weight.append(m.Intermediate(m.sum([link[b][u_] * weights[u_] for u_ in range(num_users)])))

    # Compute fairness for each user
    for u in range(num_users):
        total_gain_u.append(m.Intermediate(cumulative_gain[u]))  # Use GEKKO's Intermediate for cumulative gain

        # Loop over machines and calculate gain for user u
        for b in range(num_machines):
            # GEKKO constraint-compatible check if user u is connected to machine b
            link_value = link[b][u]

            # Use GEKKO's if3 function to avoid direct Python conditional checks
            resources_allocated_to_u_b = m.if3(total_machine_weight[b] > 0, (weights[u] / total_machine_weight[b]) * num_resources, 0)

            # Add the gain from machine b to user u, conditioned on the link value
            total_gain_u[u] += link_value * resources_allocated_to_u_b * rewards[b, u]

        # Fairness is calculated based on total gain divided by demand
        fairness_u = m.Intermediate(total_gain_u[u] / demands[u]) if demands[u] > 0 else 0
        user_fairness.append(fairness_u)

    # Compute fairness for each machine using the geometric mean of users it serves
    machine_fairness = []
    for b in range(num_machines):
        connected_user_fairness = []
        for u in range(num_users):
            connected_fairness = m.Intermediate(link[b][u] * user_fairness[u])
            connected_user_fairness.append(connected_fairness)

        if connected_user_fairness:
            product = m.Intermediate(1)
            for fairness in connected_user_fairness:
                product *= fairness
            fairness_machine = m.Intermediate(product ** (1 / len(connected_user_fairness)))
            machine_fairness.append(fairness_machine)

    # Compute total system fairness using the geometric mean of the fairness of all machines
    if machine_fairness:
        product1 = m.Intermediate(1)
        for fairness in machine_fairness:
            product1 *= fairness
        total_fairness = m.Intermediate(product1 ** (1 / len(machine_fairness)))
    else:
        total_fairness = 0
      
    total_rewards = m.Var(lb=0)  # Create GEKKO variable for total rewards

    for u in range(num_users):
        # Allocate resources to user u based on its weight relative to the total weight
        resources_allocated_to_u_b = m.if3(total_machine_weight[b] > 0, (weights[u] / total_machine_weight[b]) * num_resources, 0)

        # Calculate reward for user u based on resources allocated
        reward_contribution = link[b][u] * resources_allocated_to_u_b * rewards[b, u]
        total_rewards += reward_contribution
    
    return total_fairness, total_rewards


# Test setup
num_machines = 2
num_users = 5
num_resources = 2
demands = [100, 100, 100, 100, 100]
cumulative_gain = [0, 0, 0, 0, 0]
sinr_matrix = np.random.randint(1, num_machines * num_users, size=(num_machines, num_users))
C1 = 1
C2 = 1

print('Demands:', demands)
print('Cumulative Gain:', cumulative_gain)

link, weights, m = allocation(sinr_matrix, num_machines, num_users, num_resources, C1, C2, demands, cumulative_gain)

print('Link Matrix:', link)
print('Weights:', weights)

But i am getting following error:

 ----------------------------------------------------------------
 APMonitor, Version 1.0.1

 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            2
   Constants    :            0
   Variables    :          110
   Intermediates:           28
   Connections  :           12
   Equations    :          106
   Residuals    :           78
 
 @error: Model Expression
 *** Error in syntax of function string: Invalid element: i280>0
 
Position: 20                  
 0-((((1-int_v28))*(i280>0)))-slk_18
                    ?


---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In[50], line 120
    117 print('Demands:', demands)
    118 print('Cumulative Gain:', cumulative_gain)
--> 120 link, weights, m = allocation(rewards, num_machines, num_users, num_resources, C1, C2, demands, cumulative_gain)
    122 print('Link Matrix:', link)
    123 print('Weights:', weights)

Cell In[50], line 27, in allocation(rewards, num_machines, num_users, num_resources, C1, C2, demands, cumulative_gain)
     25 m.options.SOLVER = 1  # Change solver (1=APOPT, 3=IPOPT)
     26 m.open_folder()
---> 27 m.solve(disp=True)
     29 optimized_weights = np.array([weight.value[0] for weight in weights])
     30 optimized_link = np.ndarray((num_machines, num_users))

File ~/anaconda3/lib/python3.11/site-packages/gekko/gekko.py:2140, in GEKKO.solve(self, disp, debug, GUI, **kwargs)
   2138         print("Error:", errs)
   2139     if (debug >= 1) and record_error:
-> 2140         raise Exception(apm_error)
   2142 else: #solve on APM server
   2143     def send_if_exists(extension):

Exception: @error: Model Expression
 *** Error in syntax of function string: Invalid element: i280>0
 
Position: 20                  
 0-((((1-int_v28))*(i280>0)))-slk_18

I tried looking at the .apm file but couldnt track the issue. Following is the .apm file

Model
Variables
    int_v1 = 0, <= 1, >= 0
    int_v2 = 0, <= 1, >= 0
    int_v3 = 0, <= 1, >= 0
    int_v4 = 0, <= 1, >= 0
    int_v5 = 0, <= 1, >= 0
    int_v6 = 0, <= 1, >= 0
    int_v7 = 0, <= 1, >= 0
    int_v8 = 0, <= 1, >= 0
    int_v9 = 0, <= 1, >= 0
    int_v10 = 0, <= 1, >= 0
    v11 = 0, <= 1, >= 0.1
    v12 = 0, <= 1, >= 0.1
    v13 = 0, <= 1, >= 0.1
    v14 = 0, <= 1, >= 0.1
    v15 = 0, <= 1, >= 0.1
    v16 = 0
    v17 = 0
    v18 = 0
    v19 = 0
    v20 = 0
    v21 = 0
    v22 = 0
    v23 = 0
    v24 = 0
    v25 = 0
    v26 = 0
    v27 = 0
    int_v28 = 0.01, <= 1, >= 0
    v29 = 0
    int_v30 = 0.01, <= 1, >= 0
    v31 = 0
    int_v32 = 0.01, <= 1, >= 0
    v33 = 0
    int_v34 = 0.01, <= 1, >= 0
    v35 = 0
    int_v36 = 0.01, <= 1, >= 0
    v37 = 0
    int_v38 = 0.01, <= 1, >= 0
    v39 = 0
    int_v40 = 0.01, <= 1, >= 0
    v41 = 0
    int_v42 = 0.01, <= 1, >= 0
    v43 = 0
    int_v44 = 0.01, <= 1, >= 0
    v45 = 0
    int_v46 = 0.01, <= 1, >= 0
    v47 = 0
    v48 = 0, >= 0
    int_v49 = 0.01, <= 1, >= 0
    v50 = 0
    int_v51 = 0.01, <= 1, >= 0
    v52 = 0
    int_v53 = 0.01, <= 1, >= 0
    v54 = 0
    int_v55 = 0.01, <= 1, >= 0
    v56 = 0
    int_v57 = 0.01, <= 1, >= 0
    v58 = 0
End Variables
Intermediates
    i280=v21
    i281=v27
    i282=0
    i283=((((i282+((((int_v1)*(v29)))*(8)))+((((int_v6)*(v31)))*(2))))/(100))
    i284=0
    i285=((((i284+((((int_v2)*(v33)))*(6)))+((((int_v7)*(v35)))*(1))))/(100))
    i286=0
    i287=((((i286+((((int_v3)*(v37)))*(6)))+((((int_v8)*(v39)))*(1))))/(100))
    i288=0
    i289=((((i288+((((int_v4)*(v41)))*(6)))+((((int_v9)*(v43)))*(7))))/(100))
    i290=0
    i291=((((i290+((((int_v5)*(v45)))*(4)))+((((int_v10)*(v47)))*(9))))/(100))
    i292=((int_v1)*(i283))
    i293=((int_v2)*(i285))
    i294=((int_v3)*(i287))
    i295=((int_v4)*(i289))
    i296=((int_v5)*(i291))
    i297=1
    i298=((((((((((((i297)*(i292)))*(i293)))*(i294)))*(i295)))*(i296)))^(0.2))
    i299=((int_v6)*(i283))
    i300=((int_v7)*(i285))
    i301=((int_v8)*(i287))
    i302=((int_v9)*(i289))
    i303=((int_v10)*(i291))
    i304=1
    i305=((((((((((((i304)*(i299)))*(i300)))*(i301)))*(i302)))*(i303)))^(0.2))
    i306=1
    i307=((((((i306)*(i298)))*(i305)))^(0.5))
End Intermediates
Equations
    ((0+int_v1)+int_v6)<=1
    ((0+int_v2)+int_v7)<=1
    ((0+int_v3)+int_v8)<=1
    ((0+int_v4)+int_v9)<=1
    ((0+int_v5)+int_v10)<=1
    (((((0+((int_v1)*(v11)))+((int_v2)*(v12)))+((int_v3)*(v13)))+((int_v4)*(v14)))+((int_v5)*(v15)))<=1
    (((((0+((int_v6)*(v11)))+((int_v7)*(v12)))+((int_v8)*(v13)))+((int_v9)*(v14)))+((int_v10)*(v15)))<=1
    v16=((int_v1)*(v11))
    v17=((int_v2)*(v12))
    v18=((int_v3)*(v13))
    v19=((int_v4)*(v14))
    v20=((int_v5)*(v15))
    v22=((int_v6)*(v11))
    v23=((int_v7)*(v12))
    v24=((int_v8)*(v13))
    v25=((int_v9)*(v14))
    v26=((int_v10)*(v15))
    (((1-int_v28))*(i280>0))<=0
    ((int_v28)*(i280>0))>=0
    v29=((((1-int_v28))*(((((v11)/(i280)))*(2))))+((int_v28)*(0)))
    (((1-int_v30))*(i281>0))<=0
    ((int_v30)*(i281>0))>=0
    v31=((((1-int_v30))*(((((v11)/(i281)))*(2))))+((int_v30)*(0)))
    (((1-int_v32))*(i280>0))<=0
    ((int_v32)*(i280>0))>=0
    v33=((((1-int_v32))*(((((v12)/(i280)))*(2))))+((int_v32)*(0)))
    (((1-int_v34))*(i281>0))<=0
    ((int_v34)*(i281>0))>=0
    v35=((((1-int_v34))*(((((v12)/(i281)))*(2))))+((int_v34)*(0)))
    (((1-int_v36))*(i280>0))<=0
    ((int_v36)*(i280>0))>=0
    v37=((((1-int_v36))*(((((v13)/(i280)))*(2))))+((int_v36)*(0)))
    (((1-int_v38))*(i281>0))<=0
    ((int_v38)*(i281>0))>=0
    v39=((((1-int_v38))*(((((v13)/(i281)))*(2))))+((int_v38)*(0)))
    (((1-int_v40))*(i280>0))<=0
    ((int_v40)*(i280>0))>=0
    v41=((((1-int_v40))*(((((v14)/(i280)))*(2))))+((int_v40)*(0)))
    (((1-int_v42))*(i281>0))<=0
    ((int_v42)*(i281>0))>=0
    v43=((((1-int_v42))*(((((v14)/(i281)))*(2))))+((int_v42)*(0)))
    (((1-int_v44))*(i280>0))<=0
    ((int_v44)*(i280>0))>=0
    v45=((((1-int_v44))*(((((v15)/(i280)))*(2))))+((int_v44)*(0)))
    (((1-int_v46))*(i281>0))<=0
    ((int_v46)*(i281>0))>=0
    v47=((((1-int_v46))*(((((v15)/(i281)))*(2))))+((int_v46)*(0)))
    (((1-int_v49))*(i281>0))<=0
    ((int_v49)*(i281>0))>=0
    v50=((((1-int_v49))*(((((v11)/(i281)))*(2))))+((int_v49)*(0)))
    (((1-int_v51))*(i281>0))<=0
    ((int_v51)*(i281>0))>=0
    v52=((((1-int_v51))*(((((v12)/(i281)))*(2))))+((int_v51)*(0)))
    (((1-int_v53))*(i281>0))<=0
    ((int_v53)*(i281>0))>=0
    v54=((((1-int_v53))*(((((v13)/(i281)))*(2))))+((int_v53)*(0)))
    (((1-int_v55))*(i281>0))<=0
    ((int_v55)*(i281>0))>=0
    v56=((((1-int_v55))*(((((v14)/(i281)))*(2))))+((int_v55)*(0)))
    (((1-int_v57))*(i281>0))<=0
    ((int_v57)*(i281>0))>=0
    v58=((((1-int_v57))*(((((v15)/(i281)))*(2))))+((int_v57)*(0)))
    maximize (((1)*((((((v48+((((int_v6)*(v50)))*(2)))+((((int_v7)*(v52)))*(1)))+((((int_v8)*(v54)))*(1)))+((((int_v9)*(v56)))*(7)))+((((int_v10)*(v58)))*(9)))))+((1)*(i307)))
End Equations
Connections
    v16 = sum_1.x[1]
    v17 = sum_1.x[2]
    v18 = sum_1.x[3]
    v19 = sum_1.x[4]
    v20 = sum_1.x[5]
    v21 = sum_1.y
    v22 = sum_2.x[1]
    v23 = sum_2.x[2]
    v24 = sum_2.x[3]
    v25 = sum_2.x[4]
    v26 = sum_2.x[5]
    v27 = sum_2.y
End Connections
Objects
    sum_1 = sum(5)
    sum_2 = sum(5)
End Objects

End Model

Could you please help me in the direction where the error could be.

Improve this question

1 Answer 1

Reset to default
0

The ">0" is not needed when defining

#m.if3(condition>0,result1,result2)
m.if3(condition,result1, result2)

Here is a corrected script without the >0 in the m.if3() functions.

import numpy as np
from gekko import GEKKO

def allocation(rewards, num_machines, num_users, num_resources, C1, C2, demands, cumulative_gain):
    m = GEKKO(remote=False)
    # Decision variables
    link = m.Array(m.Var, (num_machines, num_users), lb=0, ub=1, integer=True)
    weights = m.Array(m.Var, num_users, lb=0.1, ub=1)
    # Constraints
    for u in range(num_users):
        m.Equation(m.sum(link[:, u]) <= 1)  # each user is assigned to only 1 machine
    for b in range(num_machines):
        m.Equation(m.sum(link[b, :] * weights) <= 1)  # sum of weights of users allocated to each machine should be 1
    # Compute Fairness
    fairness, reward = compute_fairness(link, weights, rewards, demands, cumulative_gain, num_resources, m)
    # Objective
    m.Maximize(C1 * reward + C2 * fairness)
    m.options.SOLVER = 1  # Change solver (1=APOPT, 3=IPOPT)
    m.open_folder()
    m.solve(disp=True)
    optimized_weights = np.array([weight.value[0] for weight in weights])
    optimized_link = np.ndarray((num_machines, num_users))
    for i in range(link.shape[0]):
        for j in range(link.shape[1]):
            optimized_link[i][j] = link[i][j][0]
    return optimized_link, optimized_weights, m

def compute_fairness(link, weights, rewards, demands, cumulative_gain, num_resources, m):
    """
    Computes the fairness of the system based on the resources allocated to users using GEKKO-compatible logic.
    """
    num_machines = link.shape[0]
    num_users = link.shape[1]
    
    total_gain_u = []
    user_fairness = []
    total_machine_weight = []
    
    for b in range(num_machines):
        total_machine_weight.append(m.Intermediate(m.sum([link[b][u_] * weights[u_] for u_ in range(num_users)])))
    # Compute fairness for each user
    for u in range(num_users):
        total_gain_u.append(m.Intermediate(cumulative_gain[u]))  # Use GEKKO's Intermediate for cumulative gain
        # Loop over machines and calculate gain for user u
        for b in range(num_machines):
            # GEKKO constraint-compatible check if user u is connected to machine b
            link_value = link[b][u]
            # Use GEKKO's if3 function to avoid direct Python conditional checks
            resources_allocated_to_u_b = m.if3(total_machine_weight[b], (weights[u] / total_machine_weight[b]) * num_resources, 0)
            # Add the gain from machine b to user u, conditioned on the link value
            total_gain_u[u] += link_value * resources_allocated_to_u_b * rewards[b, u]
        # Fairness is calculated based on total gain divided by demand
        fairness_u = m.Intermediate(total_gain_u[u] / demands[u]) if demands[u] > 0 else 0
        user_fairness.append(fairness_u)
    # Compute fairness for each machine using the geometric mean of users it serves
    machine_fairness = []
    for b in range(num_machines):
        connected_user_fairness = []
        for u in range(num_users):
            connected_fairness = m.Intermediate(link[b][u] * user_fairness[u])
            connected_user_fairness.append(connected_fairness)
        if connected_user_fairness:
            product = m.Intermediate(1)
            for fairness in connected_user_fairness:
                product *= fairness
            fairness_machine = m.Intermediate(product ** (1 / len(connected_user_fairness)))
            machine_fairness.append(fairness_machine)
    # Compute total system fairness using the geometric mean of the fairness of all machines
    if machine_fairness:
        product1 = m.Intermediate(1)
        for fairness in machine_fairness:
            product1 *= fairness
        total_fairness = m.Intermediate(product1 ** (1 / len(machine_fairness)))
    else:
        total_fairness = 0
      
    total_rewards = m.Var(lb=0)  # Create GEKKO variable for total rewards
    for u in range(num_users):
        # Allocate resources to user u based on its weight relative to the total weight
        resources_allocated_to_u_b = m.if3(total_machine_weight[b], (weights[u] / total_machine_weight[b]) * num_resources, 0)
        # Calculate reward for user u based on resources allocated
        reward_contribution = link[b][u] * resources_allocated_to_u_b * rewards[b, u]
        total_rewards += reward_contribution
    return total_fairness, total_rewards

# Test setup
num_machines = 2
num_users = 5
num_resources = 2
demands = [100, 100, 100, 100, 100]
cumulative_gain = [0, 0, 0, 0, 0]
sinr_matrix = np.random.randint(1, num_machines * num_users, size=(num_machines, num_users))
C1 = 1
C2 = 1
print('Demands:', demands)
print('Cumulative Gain:', cumulative_gain)
link, weights, m = allocation(sinr_matrix, num_machines, num_users, num_resources, C1, C2, demands, cumulative_gain)
print('Link Matrix:', link)
print('Weights:', weights)
Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.