imsl.linalg.lin_lsq_lin_constraints

lin_lsq_lin_constraints(a, b, c, bl, bu, con_type, xlb=None, xub=None, abs_tol=None, rel_tol=None, max_iter=None)

Solve a linear least-squares problem with linear constraints.

Parameters:
  • a ((M, N) array_like) – Array containing the coefficients of the M least-squares equations.
  • b ((M,) array_like) – Array containing the right-hand sides of the least-squares equations.
  • c ((K, N) array_like) – Array containing the coefficients of the K general constraints. If the problem contains no general constraints, set c = None.
  • bl ((K,) array_like) – Array containing the lower limits of the general constraints. If there is no lower limit on the i-th constraint, then bl[i] will not be referenced.
  • bu ((K,) array_like) – Array containing the upper limits of the general constraints. If there is no upper limit on the i-th constraint, then bl[i] will not be referenced.
  • con_type ((K,) array_like) – Array indicating the type of the general constraints, where con_type[i] = 0, 1, 2, 3 indicates =, <=, >= and range constraints, respectively.
  • xlb ((N,) array_like, optional) –

    Array containing the lower bounds on the variables. If there is no lower bound on the i-th variable, then xlb[i] should be set to constant UNBOUNDED_ABOVE in module imsl.constants.

    Default: The variables have no lower bounds.

  • xub ((N,) array_like, optional) –

    Array containing the upper bounds on the variables. If there is no upper bound on the i-th variable, then xlb[i] should be set to constant UNBOUNDED_BELOW in module imsl.constants.

    Default: The variables have no upper bounds.

  • abs_tol (float, optional) –

    Absolute rank determination tolerance to be used.

    Default: abs_tol = math.sqrt(numpy.finfo(numpy.float64).eps).

  • rel_tol (float, optional) –

    Relative rank determination tolerance to be used.

    Default: rel_tol = math.sqrt(numpy.finfo(numpy.float64).eps).

  • max_iter (int, optional) –

    The maximum number of add/drop iterations.

    Default: max_iter = 5 * max(M,N).
Returns:

  • A named tuple with the following fields
  • solution ((N,) ndarray) – The approximate solution of the least-squares problem.
  • residual ((M,) ndarray) – The residuals b-Ax of the least-squares equations at the approximate solution.

Notes

Function lin_lsq_lin_constraints solves linear least-squares problems with linear constraints. These are systems of least-squares equations of the form

\[Ax \cong b\]

subject to

\[\begin{split}b_l \le Cx \le b_u \\ x_l \le x \le x_u\end{split}\]

Here A is the coefficient matrix of the least-squares equations, b is the right-hand side, and C is the coefficient matrix of the constraints. The vectors \(b_l, b_u, x_l\) and \(x_u\) are the lower and upper bounds on the constraints and the variables, respectively. The system is solved by defining dependent variables \(y \equiv Cx\) and then solving the least-squares system with the lower and upper bounds on x and y. The equation \(Cx-y=0\) is a set of equality constraints. These constraints are realized by heavy weighting, i.e., a penalty method ([1]).

References

[1]Hanson, Richard J. (1986), Linear Least Squares with Bounds and Linear Constraints, SIAM J. Sci. and Stat. Computing, 7(3), 826-834.

Examples

The following problem is solved in the least-squares sense:

\[\begin{split}3x_1+2x_2+x_3=3.3 \\ 4x_1+2x_2+x_3=2.3 \\ 2x_1+2x_2+x_3=1.3 \\ x_1+x_2+x_3=1.0\end{split}\]

subject to

\[\begin{split}x_1=x_2+x_3 \le 1 \\ 0 \le x_1 \le 0.5 \\ 0 \le x_2 \le 0.5 \\ 0 \le x_3 \le 0.5\end{split}\]

The approximate solution of the least-squares problem, the residuals of the least-squares equations at the solution and the norm of the residual vector are printed.

>>> import imsl.linalg as la
>>> import numpy as np
>>> a = [[3.0, 2.0, 1.0], [4.0, 2.0, 1.0], [2.0, 2.0, 1.0],
...      [1.0, 1.0, 1.0]]
>>> b = [3.3, 2.3, 1.3, 1.0]
>>> c = [[1.0, 1.0, 1.0]]
>>> xlb = [0.0, 0.0, 0.0]
>>> xub = [0.5, 0.5, 0.5]
>>> con_type = [1]
>>> bc = [1.0]
>>> x = la.lin_lsq_lin_constraints(a, b, c, bc, bc, con_type, xlb=xlb,
...                                xub=xub)
>>> print(x.solution) 
[0.5 0.3 0.2]
>>> print(x.residual) 
[-1.   0.5  0.5  0. ]
>>> print("{0:8.6f}".format(np.linalg.norm(x.residual)))
1.224745