Solve underdetermined linear equation Ax + By = C where y is constrained by x

I am new to optimization and have been struggling to solve for variable x and y in linear equation Ax +By = C, while y is constrained by the solution of x. An example of the problem can be:

A = np.random.rand(20,100)
B = np.random.rand(20,200)
C = np.random.rand(20)

Solve for x and y so that Ax +By = C, with constraints that x is non-negative and -0.7*x0 < y0,y1 <0.7*x0, -0.7*x1 < y2,y3 <0.7*x1... ( -0.7x[i] < y[2i],y[2i+1]<0.7x[i] )

I would really appreciate if someone can recommend me a python package that solve this problem, or some way to transform my problem into a more conventional format that can be solved directly with libraries like Scipy.optimize

2 Answers

I solved the problem using CVXPY:

import cvxpy as cp
import numpy as np
m = 20
n = 100
A = np.random.rand(m,n)
B = np.random.rand(m,n)
C = np.random.rand(m,n)
d = np.random.rand(m)
# construct the problem.
x = cp.Variable(n)
y = cp.Variable(n)
z = cp.Variable(n)
objective = cp.Minimize(cp.sum_squares(A*x + B*y + C*z -d))
constaints = [0 <= x, x <= 1, y <= 0.7*x, y >= -0.7*x, z <= 0.7*x, z >= -0.7*x]
prob = cp.Problem(objective, constraints)
result = prob.solve()

It depends on what do you mean by solve, what you have described have multiple solutions and it is the interior of a polyhedral.

It is a linear programming problem if you are willing to convert the problem to say

-0.7x_0 <=y_0 <= 0.7x_0

Suppose not, consider introducing a small positive number m to make a linear program.

-0.7x_0 + m <=y_0 <= 0.7x_0 - m

You can then use scipy.optimize.linprog to solve a linear program.

If you are just interested in a particular solution, you can just set c in the objective function to be zero.

Edit:

import numpy as np
from scipy.optimize import linprog
m_size = 20
bound = 0.7
x_size = 100
y_size = 2 * x_size
small_m = 0
A = np.random.rand(m_size, x_size)
B = np.random.rand(m_size, y_size)
C = np.random.rand(m_size)
A_eq = np.hstack((A, B))
b_eq = C
sub_block = np.array([-bound, -bound ])
sub_block.shape = (2,1)
block = np.kron(np.eye(x_size), sub_block)
upper_block = np.hstack((block, - np.eye(y_size)))
lower_block = np.hstack((block, np.eye(y_size)))
A_ub = np.vstack((upper_block, lower_block))
b_ub = -small_m * np.ones( 2 * y_size)
bounds = tuple([tuple([0, None]) for i in range(x_size)] + [tuple([None, None]) for i in range(y_size)])
c = [0 for i in range(x_size+y_size)]
res = linprog(c, A_ub = A_ub, b_ub = b_ub, A_eq = A_eq, b_eq = b_eq, bounds = bounds)
x_part = res.x[:x_size]
y_part = res.x[x_size:]