Почему CVXPy или DMCP возвращают «NoneType» при добавлении простых ограничений? - PullRequest
Новая
       10

Почему CVXPy или DMCP возвращают «NoneType» при добавлении простых ограничений?

0 голосов
/ 07 ноября 2019

Недавно я наткнулся на интересный пакет / документ под названием dmcp (дисциплинированное много выпуклое программирование). Вы можете найти ссылку на репо здесь и ссылку на статью здесь .

Я хотел использовать это, чтобы найти равновесие в модели, в которой я былучусь в своих исследованиях (проблема на самом деле выпуклая, но я изо всех сил пытался перевести ее в форму DCP, поэтому я подумал, что это было бы неплохо иметь в качестве резервной копии).

Мы доказали, что равновесия всегда существуют, и я тестировал приведенный ниже код на игрушечном примере, который я решил вручную, но по какой-то причине DMCP дает сбой после нескольких итераций, так как один из CVX variable.valueстановится NoneType (я думаю).

Просматривая код репозитория DMCP, я полагаю, что ошибка происходит где-то, когда они используют CVXPy, но я не знаю достаточно, чтобы выяснить, есть ли ошибка в DMCP или в CVXPy.

Код, который я хочу получить, имеет параметры: 

num_agents = 2
num_goods = 2
adjacency_matrix = np.eye(2)
endowment = np.eye(2)
utilities = np.abs(np.eye(2)-1)
resale = np.ones(num_agents)

Переменные: 

prices = cp.Variable((num_agents,num_goods), nonneg=True)
consumption = cp.Variable((num_agents*num_agents,num_goods), nonneg=True)
resale = cp.Variable((num_agents*num_agents,num_goods), nonneg=True)
sum_x_ij = cp.Variable((num_agents,num_goods), nonneg=True)
sum_x_ji = cp.Variable((num_agents,num_goods), nonneg=True)
sum_y_ij = cp.Variable((num_agents,num_goods), nonneg=True)
sum_y_ji = cp.Variable((num_agents,num_goods), nonneg=True)

Начальные значения: 

prices.value = np.random.uniform(low=.1,high=1.1,size=(num_agents,num_goods))
consumption.value = np.ones((num_agents*num_agents,num_goods))
resale.value = np.ones((num_agents*num_agents,num_goods))

Вспомогательная функция: 

def map_to_idx(agent,neighbor):
    return (agent*num_agents)+neighbor

Цель: 

obj = 0#cp.sum(cp.multiply(utilities[:,:-1],sum_x_ij)) #+ cp.sum(cp.sum(cp.multiply(prices,sum_y_ij),axis=1))
for agent in range(num_agents):
    if utilities[agent,-1] != 1:
        obj += cp.pnorm(cp.multiply(cp.power(utilities[agent,:-1],1.0/utilities[agent,-1]),sum_x_ij[agent,:]),utilities[agent,-1])
    else:
        obj += utilities[agent,:-1]*sum_x_ij[agent,:]
    obj += prices[agent,:]*sum_y_ij[agent,:]
    for neighbor in range(num_agents):
        obj -= prices[neighbor,:]*resale[map_to_idx(agent,neighbor),:]

Ограничения: 

constraints = [prices >= 0.1, consumption >= 0, resale >= 0, sum_x_ij >= 0, sum_x_ji >= 0, sum_y_ij >= 0, sum_y_ji >= 0]
for agent in range(num_agents):
    constraints += [sum_x_ji[agent,:] + sum_y_ji[agent,:] == endowments[agent,:] + sum_y_ij[agent,:]]
    constraints += [cp.sum(sum_y_ij[agent,:]) <= resale_vect[agent]]
    lhs_budget_const = 0
    rhs_budget_const = prices[agent,:]*endowments[agent,:] + prices[agent,:]*sum_y_ij[agent,:]
    lhs_sum_xij = 0
    lhs_sum_yij = 0
    lhs_sum_xji = 0
    lhs_sum_yji = 0
    for neighbor in range(num_agents):
        lhs_budget_const += prices[neighbor,:]*consumption[map_to_idx(agent,neighbor),:]
        rhs_budget_const -= prices[neighbor,:]*resale[map_to_idx(agent,neighbor),:]
        if adjacency_matrix[agent,neighbor] == 0:
            constraints += [consumption[map_to_idx(agent,neighbor),:] == 0]
            constraints += [resale[map_to_idx(agent,neighbor),:] == 0]
        lhs_sum_xij += consumption[map_to_idx(agent,neighbor),:]
        lhs_sum_yij += resale[map_to_idx(agent,neighbor),:]
        lhs_sum_xji += consumption[map_to_idx(neighbor,agent),:]
        lhs_sum_yji += resale[map_to_idx(neighbor,agent),:]
    constraints += [lhs_budget_const == rhs_budget_const]
    constraints += [lhs_sum_xij == sum_x_ij[agent,:]]
    constraints += [lhs_sum_xji == sum_x_ji[agent,:]]
    constraints += [lhs_sum_yij == sum_y_ij[agent,:]]
    constraints += [lhs_sum_yji == sum_y_ji[agent,:]]
    constraints += [resale[map_to_idx(agent,agent),:] == 0]
constraints += [cp.sum(prices) == 1]
constraints += [cp.sum(consumption,axis=0) == cp.sum(endowments,axis=0)]

и, наконец: 

prob = cp.Problem(cp.Maximize(obj), constraints)

print()
print("minimal sets:", dmcp.find_minimal_sets(prob))
print("problem is DMCP:", dmcp.is_dmcp(prob))
print()
result = prob.solve(method='bcd',max_iter = 100)
print(result)
print()
print("p = " + str(prices.value))
print("x = " + str(sum_x_ij.value))
print("y = " + str(sum_y_ij.value))

Я знаю, что у этой программы есть решение. Когда я закомментирую ту часть цели, в которой я умножаю цены на sum_y_ij и перепродаю, то все идет гладко. Когда это присутствует в цели, я получаю следующие выходные данные / ошибки: 

Defining variables

Setting some initial guesses

Defining objective

Defining constraints


minimal sets: [[0], [1, 2, 5]]
problem is DMCP: True

max abs slack = 11.257764306622317 mu = 0.005 original objective value = 41.569236735359254 fixed objective value = 20.693015614921645 status= optimal
max abs slack = 21.155476427136296 mu = 0.005 original objective value = 81.25842921844962 fixed objective value = 60.55171440353868 status= optimal
max abs slack = 30.996381672274627 mu = 0.0075 original objective value = 120.77145211345528 fixed objective value = 99.41459375605088 status= optimal
max abs slack = 40.85265643546646 mu = 0.0075 original objective value = 160.97795623707566 fixed objective value = 138.83652995015788 status= optimal
max abs slack = 50.630481600263245 mu = 0.01125 original objective value = 200.88955474035853 fixed objective value = 177.10923936689406 status= optimal
max abs slack = 200.65501019663174 mu = 0.01125 original objective value = 442.07817550005626 fixed objective value = 319.96444147471675 status= optimal
max abs slack = 388.4237602479294 mu = 0.016875 original objective value = 671.3475469742531 fixed objective value = 540.2441656255529 status= optimal
max abs slack = 53335.2974424018 mu = 0.016875 original objective value = 54331.04746214825 fixed objective value = 35106.94351475132 status= optimal_inaccurate
max abs slack = 215936.17724686957 mu = 0.0253125 original objective value = 219254.45219690344 fixed objective value = 40873.564352143076 status= optimal_inaccurate
max abs slack = 27527746.923883904 mu = 0.0253125 original objective value = 23405815.216764793 fixed objective value = 7009847.781686075 status= optimal_inaccurate
Traceback (most recent call last):
  File "driver.py", line 64, in <module>
    solv.find_multi_convex_equilibrium(num_agents, num_goods, adjacency_matrix, endowments, utilities, resale_vect, initial_guess=None, verbose=True)
  File "/home/gabe/Graphical_Econs_with_Resale/Jain_Convex_Program/solver.py", line 123, in find_multi_convex_equilibrium
    result = prob.solve(method='bcd',max_iter = 100)
  File "/usr/local/lib/python3.5/dist-packages/cvxpy/problems/problem.py", line 289, in solve
    return solve_func(self, *args, **kwargs)
  File "/usr/local/lib/python3.5/dist-packages/dmcp-1.0.0-py3.5.egg/dmcp/bcd.py", line 72, in bcd
    result = _bcd(prob, fix_sets, max_iter, solver, mu, rho, mu_max, ep, lambd, linearize, proximal)
  File "/usr/local/lib/python3.5/dist-packages/dmcp-1.0.0-py3.5.egg/dmcp/bcd.py", line 119, in _bcd
    max_slack = np.max([np.max(np.abs(var.value)) for var in var_slack])
  File "/usr/local/lib/python3.5/dist-packages/dmcp-1.0.0-py3.5.egg/dmcp/bcd.py", line 119, in <listcomp>
    max_slack = np.max([np.max(np.abs(var.value)) for var in var_slack])
TypeError: bad operand type for abs(): 'NoneType'

Эта ошибка даже происходит, если я даже включаю в цель единичные цены [0,0] * перепродажу [0,0].

Любой совет будет очень полезным! Спасибо!

...