0

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

hpaulj
  • 221,503
  • 14
  • 230
  • 353
Susie
  • 287
  • 1
  • 13

2 Answers2

0

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:]
Siong Thye Goh
  • 3,518
  • 10
  • 23
  • 31
  • Hi, thank you for your comments and for suggesting the Scipy function. I looked into the function carefully but still don't know how to solve my particular problem with the function. I wonder if you can be more specific with your answer and maybe show how the problem can be plugged in with the suggested function? I am interested in finding one feasible solution, and also C is non-zero, and -0.7x[i] – Susie Aug 28 '18 at 14:39
  • I have included the code, though I doubt it is understandable if you do not get what is going on. So let's figure out what is stopping you. Can you first see that it is a linear programming problem in the first place? note that there are two "C/c", the random C that you generated and also the objective function c in the optimization model. Or are you having difficulty in constructing the matrices? – Siong Thye Goh Aug 28 '18 at 16:34
  • This method doesn't appear to be very robust? Some constraints seem to be verified afterward ... – Susie Aug 28 '18 at 20:06
  • I am checking the code is working explicitly. your solution seems to be assumed that <= rather than <. I have modified the code accordingly though I think your code is cooler using cvxpy – Siong Thye Goh Aug 29 '18 at 00:43
0

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()
Susie
  • 287
  • 1
  • 13