0

I'm hoping to just get some assistance conceptually about how to solve a linear system of equations with penalty functions. Example code is at the bottom.

Let's say I'm trying trying to do a fit of this equation:

Ie=IaX+IbY+IcZ

where Ie, Ia, Ib, and Ic are constants, and X,Y,Z are variables I could easily solve this system of equations using scipy.least_squares, but I want to constrain the system with 2 contraints.

1. X+Y+Z=1
2. X,Y,Z > 0 and X,Y,Z < 1

To do this, I then modified the above function.

Ie-Ic=X(Ia-Ic)+Y(Ib-Ic) where X+Y+Z=1 I solved for Z
therefore
Ax=B where A=[Ia-Ic,Ib-Ic] and B=[Ie-Ic] given bounds (0,1)

This solves the 2nd criteria of X,Y,Z > 0 and X,Y,Z < 1, but it does not solve the 1st. To resolve the first issue, an additional constraint needs to be made, where X+Y<1, and this I don't quite know how to do.

So I presume least_squares has some built in penalty function for it's bounds. I.E.

chi^2=||A-Bx||^2+P
where P is the conditions, if X,Y,Z>0 then P = 10000
thus, giving high chi squared values and thus setting the constraints

I don't know how I can add a further condition. So if X+Y<1 then P=10000 or something similar along those lines.

In short, least_squares enables you to set bounds on the individual values, but I'd like to set some further constraints, and I don't quite know how to do this with least_squares. I've seen additional inequality constraint options in scipy.minimize, but I don't quite know how to apply that to a linear system of equations with the format Ax=B.

So as an example code, lets say I've already done the calculations and obtained my A matrix of constants and my B vector. I use least squares, get my values for X and Y, and can calculate Z since X+Y+Z=1. The issue here is in my minimization, I did not set a constraint that X+Y<1, so in some cases you can actually get values where X+Y>1. So I'd like to find a method where I can set that additional constraint, in addition to the bounds constraint for the individual variables:

Ax=np.array([[1,2],[2,4],[3,4]])
B=np.array([0,1,2])

solution=lsq_linear(Ax,B,lsq_solver='lsmr',bounds=(0,1))

X=solution.x[0]
Y=solution.x[1]
Z=1-sum(solution.x)

If minimize is the solution here, can you please show me how to set it up given the above matrix of A and array of B?

Any advice, tips, or help to point me in the right direction is greatly appreciated!

Edit: So I found something similar on here: Minimizing Least Squares with Algebraic Constraints and Bounds So I thought I'd apply it to my case, but I don't think I've been able to apply it properly.


Ax=np.array([[1,2],[2,4],[3,4]])
B=np.array([0,1,2])
def fun(x,a1,a2,y):
   fun_output=x[0]*a1+x[1]*a2
   return np.sum((fun_output-y)**2)


cons = [{"type": "eq", "fun": lambda x: x[0] + x[1] - 1}]
bnds = [(0, 1), (0, 1)] 
xinit = np.array([1, 1]) 
solution=minimize(fun,args=(Ax[:,0],Ax[:,1], B), x0=xinit, bounds=bnds, constraints=cons)
solution_2=lsq_linear(Ax,B,bounds=(0,1))
print(solution.x)
print(solution_2.x)

Issue is, the output of this differs from lsq_linear, and I almost always get a very close to zero value for Z regardless of what the input arrays are. I don't know if I'm setting this up/understanding this correctly.

samman
  • 559
  • 10
  • 29
  • If I understand you correctly, you basically just want to fit a function subject to some additional constraints on the (unknown) function coefficients (= the variables to optimize). If so, then formulating this as a constrained optimization problem and solving it by means of scipy.optimize.minimize in the same vein as your linked answer is probably the best approach. That being said, could you please provide your given constants Ie, Ia etc. and the data points against which you want to fit your equation? – joni Nov 15 '22 at 10:15
  • Correct. My example code has exactly what you asked. I have my linear function Ax=B where A is the constants, x is the variables, and B is the solution. My issue now is that my minimize function, provides a different output from my lsq_linear. In fact, the sum of the solutions in the minimize route gives me a value above 1, which should not even be possible given my constraint. This is why I suspect the minimize route I have setup is correct, but I think the way I have written the code might be problematic. – samman Nov 15 '22 at 14:56

1 Answers1

1

Your initial guess xinit is not feasible and doesn't satisfy your constraint.

IMO, solving the initial problem directly as a constrained nonlinear optimization problem (NLP) instead of rewriting it is the easier approach. Assuming you have all the data points Ia, Ib, Ic and Ie (you didn't provide all of them), you can use the following code snippet which is in the same vein as the linked answer of mine in your question.

from scipy.optimize import minimize
import numpy as np

def fun_to_fit(coeffs, *args):
    x, y, z = coeffs
    Ia, Ib, Ic = args
    return Ia*x + Ib*y + Ic*z

def objective(coeffs, *args):
    Ia, Ib, Ic, Ie = args
    residual = Ie - fun_to_fit(coeffs, Ia, Ib, Ic)
    return np.sum(residual**2)

# Constraint: x + y + z == 1
con = [{'type': 'eq', 'fun': lambda coeffs: np.sum(coeffs) - 1}]

# bounds
bounds = [(0, 1), (0, 1), (0, 1)]

# initial guess (fulfils the constraint and lies within the bounds)
x0 = np.array([0.25, 0.5, 0.25])

# your given data points
#Ia = np.array(...)
#Ib = np.array(...)
#Ic = np.array(...)
#Ie = np.array(...)

# solve the NLP
res = minimize(lambda coeffs: objective(coeffs, Ia, Ib, Ic, Ie), x0=x0, bounds=bounds, constraint=con)
joni
  • 6,840
  • 2
  • 13
  • 20
  • Ah so it's just the original data setup that's an issue, not the way my code is organized (those were just random numbers I made for example code). I was able to resolve the issue finally by using linear contraint. ```linear_constraint = LinearConstraint([1, 1], [0], [1])``` and for minimize using ```method='trust-constr'```. These seem like simple quick solutions over converting the system to a NLS though. – samman Nov 15 '22 at 18:37
  • @samman If my answer was helpful, I'd appreciate it if you'd accept and/or upvote it, though. – joni Nov 16 '22 at 09:28
  • I apologize for the late reply. Been a busy few days. Thank you for your answer, it was helpful! – samman Nov 21 '22 at 16:36