Search code examples
constraintsnonlinear-optimizationmixed-integer-programminggams-math

How to limit the count of nonzero variables in GAMS mixed-integer nonlinear programming?


The background of the problem:

There are 25 candidates on Neo blockchain that receive votes. Every voter except me has voted. The candidates ranked 1st to 7th each will give 200 dollars, proportionally to people who vote for them. And the 8th to 21st will each distribute 100 dollars. For example, I cast 100 votes to the currently 8th candidate, pushing him/her to 7th position, and other people voted 900 votes for him, so I get 200*100/(900+100)=20 dollars from the candidate.

I know the voting result yielded by all the other voters, and I can additionally cast N votes to n candidates at maximum, because I have only n agents as my delegates. That is, my voting vector should have no more than n nonzero values, and the sum of my voting vector should equal to N. How should I vote for highest profit?

("Dollar" is written as "GAS" in the following codes. And the number of votes I can cast is written as "NEO" in the codes)

My GAMS codes:

These codes makes an error: Endogenous $ operation not allowed
(I have not figured out how to deal with the change of position because of my vote. This may be very difficult. ) (LINDO solver gives the best results for now)

Set i "candidate index" /
    0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24
    /;

Parameters
    V(i) "other people's votes" /
        0 11409980
        1 3024833
        2 2200734
        3 2040928
        4 2032100
        5 2017067
        6 2010927
        7 2007339
        8 2001625
        9 1008817
        10 1000025
        11 1000000
        12 641913
        13 396926
        14 392171
        15 390364
        16 383919
        17 383722
        18 382447
        19 382310
        20 381850
        21 923
        22 272
        23 180
        24 0
    /
    NEO  "our NEOs" /1162/
    num_agents  "our num of agents" /2/
    reward_factor(i) "how much gas is given by the candidate at each rank" /
        0 200
        1 200
        2 200
        3 200
        4 200
        5 200
        6 200
        7 100
        8 100
        9 100
        10 100
        11 100
        12 100
        13 100
        14 100
        15 100
        16 100
        17 100
        18 100
        19 100
        20 100
        21 0
        22 0
        23 0
        24 0
    /;
    
Integer Variable vote(i) "voting vector that needs to be optimized";
vote.up(i) = NEO;
Free Variable GAS_reward "total GAS reward";

Binary Variable is_nonzero(i) "whether an element in the voting vector is zero";
Integer Variable count_nonzero "count of nonzero values in vote(i)";


Equations
    GAS_reward_eqn "the reward according to our voting vector v"
    sum_v_eqn "sum of voting vector equals num of our NEOs"
    is_nonzero_eqn(i) "whether the element vote(i) is not zero"
    num_agents_eqn "sum of nonzero element in voting vector does not exceed num of our agents";

GAS_reward_eqn..  GAS_reward =e= sum(i, reward_factor(i)*vote(i)/(V(i)+vote(i)+0.000001));
sum_v_eqn..           NEO =e= sum(i, vote(i));
is_nonzero_eqn(i)..  is_nonzero(i) =e= 1 $ (vote(i));
num_agents_eqn..  sum(i, is_nonzero(i)) =l= num_agents;

Model NEOBurger /all/;
Option MINLP=lindo;
Solve NEOBurger using MINLP maximizing GAS_reward;

(where the code is_nonzero_eqn(i).. is_nonzero(i) =e= 1 $ (vote(i)); is nonlinear!)

If you need extra test cases and references, here are the Python (3.7 or 3.8) codes for the problem strategy_multi_agents.py

# Good profit
# Can somehow deal with sorting problem
# Cannot restrict num of agents
# Cannot restrict to integer
from typing import Union, List

import numpy as np
np.set_printoptions(suppress=True)
from scipy.optimize import minimize

class BurgerNEOMultiAgentStrategy:
    """
    To find the best distribution of our GAS assuming that we do not affect the ranking
    """
    consensus_rank, consensus_gain = 7, 200
    council_rank, council_gain = 21, 100
    
    def __init__(self, candidate_votes: List[int], our_NEOs: int, num_of_agents: int = 2):
        """
        :param candidate_votes: should be sorted in descending order. All the results are based on the descending order.
        """
        self.len_candidate_votes = len(candidate_votes)
        if self.len_candidate_votes >= self.council_rank:
            self.candidate_votes = sorted(candidate_votes, reverse=True)
            if self.candidate_votes != candidate_votes:
                raise ValueError('You did not input a descending list of candidate_votes.')
            self.candidate_votes = np.array(self.candidate_votes, dtype=np.int_)  # descending order
        else:
            raise ValueError(f'len(candidate_votes) should be >= {self.council_rank}. Got {self.len_candidate_votes}.')
        self.our_NEOs = np.int_(our_NEOs)
        reward_factor = np.zeros(self.len_candidate_votes, dtype=np.int_)
        reward_factor[0:self.consensus_rank] = self.consensus_gain
        reward_factor[self.consensus_rank:self.council_rank] = self.council_gain
        self.reward_factor = reward_factor
        self.num_of_agents = num_of_agents
    
    def calc_GAS_gain(self, voting_vector: np.ndarray) -> Union[np.float_, np.int_]:
        new_voting_result = voting_vector + self.candidate_votes
        argsort = np.flip(new_voting_result.argsort())
        voting_vector = voting_vector[argsort]
        new_voting_result = new_voting_result[argsort]
        
        return np.sum(self.reward_factor * voting_vector / (new_voting_result + 1e-5))

    def calc_voting_vector(self):
        def sum_is_our_NEO(voting_vector: np.ndarray):
            return np.sum(voting_vector) - self.our_NEOs
        starting_point = self.our_NEOs / self.len_candidate_votes * np.ones(self.len_candidate_votes)
        constraints = [
            {'type': 'ineq', 'fun': lambda voting_vector: voting_vector},  # all greater than or equals 0
            {'type': 'eq','fun': sum_is_our_NEO},  # sum is our_NEO
            # {'type': 'ineq', 'fun': lambda x: -np.abs(x - np.round(x))},  # all integers
            # {'type': 'ineq', 'fun':
            #     lambda x: - np.count_nonzero(np.invert(np.isclose(x, np.zeros(self.len_candidate_votes), atol=1e-3/self.len_candidate_votes))) + self.num_of_agents},  # num of agents
        ]
        result = minimize(lambda voting_vector: (-self.calc_GAS_gain(voting_vector)),
                          starting_point, constraints=constraints, tol=1e-15, options={'maxiter': 1e9})
        # return np.round(result.x).astype(np.int_), -result.fun, np.sum(np.round(result.x).astype(np.int_))
        return result.x, -result.fun, np.sum(result.x)

if __name__ == '__main__':
    print(
        BurgerNEOMultiAgentStrategy(
        [11409980, 3024833, 2200734, 2040928, 2032100, 2017067, 2010927, 2007339, 2001625, 1008817, 1000025, 1000000, 641913, 396926, 392171, 390364, 383919, 383722, 382447, 382310, 381850, 923, 272, 180, 0],
        1162, 2).calc_voting_vector()
    )
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,99, 75,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 0,0,0,0, 0,0], 50, 2).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,99, 75,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 0,0], 25, 2).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,99, 75,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 1,1], 25, 2).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,5, 5,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 1,1], 25, 2).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,5, 5,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 1,0], 25, 2).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,5, 5,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 1,0], 1000, 2).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,5, 5,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 0,0], 26, 2).calc_voting_vector())
    print()

# Restricted to integer
# Restricted num of agents
# Not robust (sometimes constraints are broken)
# Not good profit
# Cannot deal with sorting
from typing import List
import random
from gekko import GEKKO
from strategy_multi_agents import BurgerNEOMultiAgentStrategy as ScipyStrategy

class BurgerNEOMultiAgentStrategy:
    """
    To find the best distribution of our GAS assuming that we do not affect the ranking
    """
    consensus_rank, consensus_gain = 7, 200
    council_rank, council_gain = 21, 100
    
    def __init__(self, candidate_votes: List[int], our_NEOs: int, num_of_agents: int = 2):
        """
        :param candidate_votes: should be sorted in descending order. All the results are based on the descending order.
        """
        self.model = GEKKO(remote=False, name='neoburger-strategy')
        self.model.options.IMODE = 3
        self.model.options.MAX_ITER = 10000
        self.model.options.COLDSTART = 1
        self.model.options.AUTO_COLD = 1

        self.len_candidate_votes = len(candidate_votes)
        if self.len_candidate_votes >= self.council_rank:
            self.candidate_votes = sorted(candidate_votes, reverse=True)
            if self.candidate_votes != candidate_votes:
                raise ValueError('You did not input a descending list of candidate_votes.')
            self.candidate_votes = [self.model.Param(can) for can in self.candidate_votes]  # descending order
        else:
            raise ValueError(f'len(candidate_votes) should be >= {self.council_rank}. Got {self.len_candidate_votes}.')
        self.our_NEOs = self.model.Param(our_NEOs, name='our_NEOs')
        reward_factor = [self.consensus_gain] * self.consensus_rank \
                      + [self.council_gain] * (self.council_rank - self.consensus_rank) \
                      + [0] * (self.len_candidate_votes - self.council_rank)
        self.reward_factor = [self.model.Param(r) for r in reward_factor]
        self.num_of_agents = self.model.Param(num_of_agents, name='num_of_agents')
        
        # voting_vector = [self.our_NEOs / self.len_candidate_votes] * self.len_candidate_votes
        voting_vector = ScipyStrategy(candidate_votes, our_NEOs, num_of_agents).calc_voting_vector()[0]
        voting_vector = [max(0, round(v)) for v in voting_vector]
        self.voting_vector = [self.model.Var(v, lb=0, ub=our_NEOs, integer=True) for v in voting_vector]
        
    def calc_GAS_gain(self):
        # new_vote_result = [v+V for v, V in zip(self.voting_vector, self.candidate_votes)]
        # sorted_vote_result = [(v, V) for v, V in zip(self.voting_vector, new_vote_result)]  # sorted(, key=lambda x: x[1].value)
        #
        # sorted_voting_vector = [v for v, V in sorted_vote_result]
        # sorted_vote_result = [V for v, V in sorted_vote_result]
        # return sum([v*r/(V+v+0.0001) for v, V, r in zip(sorted_voting_vector, sorted_vote_result, self.reward_factor)])
        return sum([v*r/(V+v+1) for v, V, r in zip(self.voting_vector, self.candidate_votes, self.reward_factor)])

    def calc_voting_vector(self):
        # constraints
        # sum(voting_vector) == our_NEOs
        self.model.Equation(sum(self.voting_vector) == self.our_NEOs)
        # agent number
        nonzero_vote = sum([self.model.if3(v - 0.1, 0, 1) for v in self.voting_vector])
        # self.model.Equation(nonzero_vote >= 0)
        self.model.Equation(nonzero_vote <= self.num_of_agents)
        
        # objective
        self.model.Maximize(self.calc_GAS_gain())
        
        # self.model.options.RTOL = 1e-14
        self.model.options.OTOL = 1e-14
        self.model.options.REQCTRLMODE = 1
        self.model.solve(disp=False, debug=False, GUI=False)
        return self.voting_vector, -self.model.options.OBJFCNVAL


if __name__ == '__main__':
    print(
        BurgerNEOMultiAgentStrategy(
        [11409980, 3024833, 2200734, 2040928, 2032100, 2017067, 2010927, 2007339, 2001625, 1008817, 1000025, 1000000, 641913, 396926, 392171, 390364, 383919, 383722, 382447, 382310, 381850, 923, 272, 180, 0],
        1162, 6).calc_voting_vector()
    )
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,99, 75,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 0,0,0,0, 0,0], 50, 6).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,99, 75,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 0,0], 25, 6).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,99, 75,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 1,1], 25, 6).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,5, 5,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 1,1], 25, 6).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,5, 5,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 1,0], 25, 6).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,5, 5,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 1,0], 1000, 6).calc_voting_vector())
    print(BurgerNEOMultiAgentStrategy([100,100,100,100,100, 100,5, 5,5,5, 5,5,5,5,5, 5,5,5,5,5, 5, 1,1,1,1, 0,0], 26, 6).calc_voting_vector())
    print()

And I also have codes (not shown here) which makes everything correct for only 1 or 2 agents. This is achieved not by optimization, but by manual mathematical derivation.


Solution

  • I would write this as:

    * This is not using proper GAMS syntax:
    *is_nonzero_eqn(i)..  is_nonzero(i) =e= 1 $ (vote(i));
    * a linear version is:
    is_nonzero_eqn(i)..   vote(i) =l= is_nonzero(i)*vote.up(i);
    

    This is a standard MIP formulation for this. Using:

    Model NEOBurger /all/;
    Option MINLP=baron, optcr = 0;
    Solve NEOBurger using MINLP maximizing GAS_reward;
    

    we see:

    GAMS/BARON       36.1.0 r2c0a44a Released Aug 02, 2021 WEI x86 64bit/MS Window
    
    ===========================================================================
     BARON version 21.1.13. Built: WIN-64 Wed Jan 13 16:04:12 EST 2021
    
     BARON is a product of The Optimization Firm.
     For information on BARON, see https://minlp.com/about-baron
    
     If you use this software, please cite publications from
     https://minlp.com/baron-publications, such as: 
    
     Kilinc, M. and N. V. Sahinidis, Exploiting integrality in the global
     optimization of mixed-integer nonlinear programming problems in BARON,
     Optimization Methods and Software, 33, 540-562, 2018.
    ===========================================================================
     This BARON run may utilize the following subsolver(s)
     For LP/MIP/QP: ILOG CPLEX                                      
     For NLP: MINOS, SNOPT, GAMS external NLP, IPOPT, FILTERSD, FILTERSQP
    ===========================================================================
     Doing local search
     Solving bounding LP
     Preprocessing found feasible solution with value 0.203660597686E-001
     Starting multi-start local search
     Done with local search
    ===========================================================================
      Iteration    Open nodes         Time (s)    Lower bound      Upper bound
    *         1             1             0.50     0.303390         3.71066    
              1             1             0.56     0.303390         0.303737    
    *         2             2             0.62     0.303405         0.303737    
    *         2             2             0.69     0.303428         0.303737    
    *         4             2             0.83     0.303473         0.303737    
    *         4             1             0.88     0.303639         0.303737    
    *         5             1             0.92     0.303681         0.303737    
              5             0             0.92     0.303681         0.303681    
    
     Calculating duals
    
                             *** Normal completion ***            
    
     Wall clock time:                     4.05
     Total CPU time used:                 0.94
    
     Total no. of BaR iterations:       5
     Best solution found at node:       5
     Max. no. of nodes in memory:       2
     
     All done
    ===========================================================================
    
    Solution      = 0.30368112407161  found at node 5
    Best possible = 0.303681125072
    Absolute gap  = 1.00039027062238E-9  optca = 1E-9
    Relative gap  = 3.29421287011233E-9  optcr = 1E-9
      (Note that BARON uses a different formula to compute the relative gap as
       was used for the above reported value.)
    

    If we display the pertinent variables we see:

    ----     90 VARIABLE vote.L  voting vector that needs to be optimized
    
    19 466.000,    20 696.000
    
    
    ----     90 VARIABLE is_nonzero.L  whether an element in the voting vector is zero
    
    19 1.000,    20 1.000
    

    This is correct behavior.

    Note that MINLP = NLP + MIP so we can (and should) use all the tools and tricks we know about MIP modeling when writing MINLP models.