Trying to use scipy to solve system of linear equations but having trouble

Question

Trying to solve for Nash Equilibrium

Array I am trying to solve is:

[[0, 1, 1, 1, 1],
[-1, 0, 1, 1, 1],
[-1, -1, 0, 1, 1],
[-1, -1, -1, 0, 1],
[-1, -1, -1, -1, 0]]

I know there is an equilibrium at [0][0]

What I want the output to produce is an array [1, 0, 0, 0, 0] showing the optimal strategy is move at index 0

Looking at slides, we are supposed to set up a system of linear equations that maximize V such that

0 * x1 + 1 * x2 + 1 * x3 + 1 * x4 + 1 * x5 - v >=0 for row 1, repeat for all rows. Also x1 through xn all add up to 1.

I try to set up a system of linear equations using scipy keep getting either the wrong message or an error. Right now I get "Optimization failed: The problem is infeasible. (HiGHS Status 8: model_status is Infeasible; primal_status is At lower/fixed bound"

But previously I got an array of [0 0 0 0 0] which is obviously nonsensical as a player has to make a move and also I once got an array of [1 -1 1 -1 1] which DOES solve the system of linear equations but I don't want xn to be in the negatives as these represent probabiliites

    def solve_nash_equilibrium(NashArray):
        n = NashArray.shape[0]

        # Coefficients for the objective function (maximizing v)
        c = np.zeros(n + 1)
        c[-1] = 1  # Coefficient for v (to maximize)

        # Coefficients for the equality constraints (M * x + v == 0)
        A_eq = np.hstack((NashArray, -np.ones((n, 1))))  # Coefficients for x and v
        b_eq = np.zeros(n)  # Right-hand side of the equalities

        # Additional equality constraint: sum of probabilities equals 1
        A_eq_sum = np.ones((1, n + 1))
        b_eq_sum = np.ones(1)

        # Stack the new equality constraint with the existing ones
        A_eq = np.vstack((A_eq, A_eq_sum))
        A_eq[-1][-1] = 0
        print(A_eq)
        b_eq = np.hstack((b_eq, b_eq_sum))
        print(b_eq)

        # Bounds for each variable (x and v)
        x_bounds = [(0, 1)] * n  # Players can choose any probability distribution
        v_bounds = (0, 2)  # v can take any value

        # Concatenate bounds for x and v
        bounds = x_bounds + [v_bounds]

        # Solve the linear programming problem

        result = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=[*x_bounds, v_bounds])

        # Extract the solution
        if result.success:
            strategy = result.x[:-1]  # Extract mixed strategy for each player
            max_v = 1  # Maximum value of v (the value of the game)

            print("Mixed strategy for each player:", strategy)
            print("Value of the game (maximum v):", max_v)
            return strategy
        else:
            print("Optimization failed:", result.message)
            return None

Answer 1

Use milp instead of linprog; also notably, maximization requires a coefficient of -1. This probably still doesn't make much sense due to v=0, but that's likely a problem in your definitions (is v supposed to be a vector?)

import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint


def solve_nash_equilibrium(nash_array: np.ndarray) -> tuple[
    np.ndarray,
    float,
]:
    m, n = nash_array.shape

    # sum(x) = 1
    x_sum = LinearConstraint(
        A=np.concatenate((np.ones(n), (0,))),
        lb=1, ub=1,
    )
    
    # nash_array.x - v >= 0
    nash_constraint = LinearConstraint(
        A=np.hstack((nash_array, -np.ones((m, 1)))),
        lb=0,
    )

    result = milp(
        # Maximize v
        c=np.concatenate((np.zeros(n), (-1,))),
        integrality=np.zeros(n + 1, dtype=np.uint8),
        bounds=Bounds(lb=0),
        constraints=(x_sum, nash_constraint),
    )
    if not result.success:
        raise ValueError(result.message)
    x, (v,) = np.split(result.x, (-1,))
    return x, v


def demo() -> None:
    x, v = solve_nash_equilibrium(
        nash_array=np.array((
            ( 0,  1,  1,  1,  1),
            (-1,  0,  1,  1,  1),
            (-1, -1,  0,  1,  1),
            (-1, -1, -1,  0,  1),
            (-1, -1, -1, -1,  0),
        )),
    )
    print('x =', x)
    print('v =', v)


if __name__ == '__main__':
    demo()

x = [0. 0. 0. 0. 1.]
v = -0.0