Solving LP using PuLP with numpy array

Question

I want to use numpy matrix with PuLP to set constraints.

I've a 2x4x4 numpy matrix and I want to use this matrix for constraints but the problem I've is how to use this. Actually I'm facing problem in indexing as I've to loop over all variables and fix the contraints. These are the matrices.

P = np.array([[[0.7, 0.3,0,0],
               [0,0.7,0.3,0],
               [0,0,0.6,0.4],
               [0,0,0,1]],
              [[0.7,0.3,0,0],
               [0.7,0.3,0,0],
               [0.7,0.3,0,0],
               [0.7,0.3,0,0]]])
C = np.array([[100,80,50,10],[-100,-100,-100,-100]])
beta = 0.9

P matrix is probability matrix and second one is cost matrix. Every 4x4 matrix depicts the transition probability from one state to another. and my constraint is

$\forall i:\ V_i - \beta \sum_{j=i}^{N}P[d][i][j]]V_j \ge C[i][d]$

Here V is variable.

kabdulla kabdulla · Accepted Answer · 2019-06-08T14:19:00

I'm going to assume two things;

That in that last constraint you mean C[d][i] on right-hand side, rather than C[i][d]... because P.shape[0] = d = 2, and C.shape[0] = 2.
That you are wanting the constraints to be for all d, as well as for all i.

Assuming the above, the following should do what you want:

from pulp import *
import numpy as np 

P = np.array([[[0.7, 0.3,0,0],
               [0,0.7,0.3,0],
               [0,0,0.6,0.4],
               [0,0,0,1]],
              [[0.7,0.3,0,0],
               [0.7,0.3,0,0],
               [0.7,0.3,0,0],
               [0.7,0.3,0,0]]])

C = np.array([[100,80,50,10],[-100,-100,-100,-100]])

beta = 0.9

set_D = range(0, P.shape[0])
set_I = range(0, P.shape[1])

# Generate proble, & Create variables
prob = LpProblem("numpy_constraints", LpMinimize)
V = pulp.LpVariable.dicts("V", set_I, cat='Continuous')

# Make up an objective, let's say sum of V_i
prob += lpSum([V[i] for i in set_I])

# Apply constraints
for d in set_D:
    for i in set_I:
        prob += V[i] - beta*lpSum([P[d][i][j]*V[j] for j in set_I]) >= C[d][i]

# Solve problem
prob.solve()

# Print results:
V_soln = np.array([V[i].varValue for i in set_I])
print (("Status:"), LpStatus[prob.status])
print("V_soln: ")
print(V_soln)

With which I get the following. I've not checked your constraints are satisfied but they should be.

Status: Optimal
V_soln: 
[690.23142 575.50231 492.35502 490.23142]

Solving LP using PuLP with numpy array

1 Answers