# How do I constrain a fitted curve through specific points like the origin in MATLAB?

248 views (last 30 days)

Show older comments

MathWorks Support Team
on 3 May 2012

Edited: MathWorks Support Team
on 18 Oct 2023

### Accepted Answer

MathWorks Support Team
on 17 Oct 2023

Edited: MathWorks Support Team
on 18 Oct 2023

As of MATLAB R2023b, constraining a fitted curve so that it passes through specific points requires the use of a linear constraint. Neither the 'polyfit' function nor the Curve Fitting Toolbox allows specifying linear constraints. Performing this operation requires the use of the 'lsqlin' function in the Optimization Toolbox.

Consider the data created by the following commands:

c = [1 -2 1 -1];

x = linspace(-2,4);

y = c(1)*x.^3+c(2)*x.^2+c(3)*x+c(4) + randn(1,100);

plot(x,y,'.b-')

You can view the unconstrained fit to a third-order polynomial (using POLYFIT) via:

hold on

c = polyfit(x,y,3);

yhat = c(1)*x.^3+c(2)*x.^2+c(3)*x+c(4);

plot(x,yhat,'r','linewidth',2)

However, if you wish to constrain the fit to go through a specific point, for example (x0, y0) where:

x0 = 1;

y0 = 10;

use the LSQLIN function in the Optimization Toolbox to solve the linear least-squares problem with a linear constraint, as in the following example:

x = x(:); %reshape the data into a column vector

y = y(:);

% 'C' is the Vandermonde matrix for 'x'

n = 3; % Degree of polynomial to fit

V(:,n+1) = ones(length(x),1,class(x));

for j = n:-1:1

V(:,j) = x.*V(:,j+1);

end

C = V;

% 'd' is the vector of target values, 'y'.

d = y;

%%

% There are no inequality constraints in this case, i.e.,

A = [];

b = [];

%%

% We use linear equality constraints to force the curve to hit the required point. In

% this case, 'Aeq' is the Vandermoonde matrix for 'x0'

Aeq = x0.^(n:-1:0);

% and 'beq' is the value the curve should take at that point

beq = y0;

%%

p = lsqlin( C, d, A, b, Aeq, beq )

%%

% We can then use POLYVAL to evaluate the fitted curve

yhat = polyval( p, x );

%%

% Plot original data

plot(x,y,'.b-')

hold on

% Plot point to go through

plot(x0,y0,'gx','linewidth',4)

% Plot fitted data

plot(x,yhat,'r','linewidth',2)

hold off

##### 2 Comments

Sandip Kumar
on 29 Jan 2015

Edited: MathWorks Support Team
on 23 Feb 2022

K E
on 27 Mar 2015

Edited: K E
on 27 Mar 2015

Thanks so much. Saved me a ton of time. Note that in R2015a, returns an error:

Warning: The trust-region-reflective algorithm can handle bound constraints only;

using active-set algorithm instead.

> In lsqlin (line 304)

What is the suggested workaround? Or perhaps it's not necessary; the function still returns the answer (p).

### More Answers (2)

Are Mjaavatten
on 28 Nov 2015

Edited: MathWorks Support Team
on 15 May 2023

##### 2 Comments

Are Mjaavatten
on 10 Feb 2016

Hello Sheila,

Yes, you may specify more than one point, as demonstrated by example 1 in the help text. But note that if you specify a polynomial degree <= the number of contraints, the data points are ignored.

Example:

x = linspace(0,4,200)';y = sin(pi*x)+ 0.1*randn(200,1);

p = polyfix(x,y,3,[0,2,4],[0,0,0]);plot(x,y,'.',x,polyval(p,x));

You get a better result by using a polynomial of degree 7:

p = polyfix(x,y,7,[0,2,4],[0,0,0]);plot(x,y,'.',x,polyval(p,x));

However: Take care when you use polynomials of such high degree. They may not behave as you expect, especially outside the range of the fixed points.

John D'Errico
on 27 Mar 2015

Well, actually, this problem can be solve quite trivially in a way that does not require lsqlin, or anything fancy. There really is no requirement for a linear constraint, since the constraint will be that the constant parameter in the polynomial model is just zero. So use backslash to estimate the model coefficients, and just leave out the constant term. Note that this works ONLY to force a point through the origin, since the point (x,y) = (0,0) in a polynomial model implies that the constant term MUST be identically zero.

So, given a model of the form

y = a1*x^3 + a*x^2 + a3*x

we can fit it simply as

a123 = [x.^3, x.^2, x]\y;

where x and y are column vectors of data.

You can even use polyval to evaluate the polynomial, simply by appending a zero to the end of that vector f coefficients, thus

P = [a123;0];

If you wanted to force the mode to pass through the point (x0,y0), where these values are not the origin, then we need to work slightly more, but not a lot. While you can do so with lsqlin, if you don't have that tool (or my own LSE, as found on the file exchange) just do this simple transformation:

xtr = x - x0;

a123 = [xtr.^3, xtr.^2, xtr]\(y - y0);

Essentially, we have now fit the model

y - y0 = a1*(x - x0)^3 + a*(x - x0)^2 + a3*(x - x0)

You can see that when x == x0, y must be predicted as y0.

Evaluation of this model is just as easy. You can still use polyval, just remember to subtract off x0 from x first.

P = [a123;y0];

ypred = polyval(P,xpred - x0);

What happens when we have a more complicated problem, where maybe you need to force TWO points on the curve? Well, then, you can use either lsqlin or my own LSE. (Note that LSE runs in R2015a with no warning messages.)

Or, if you have a trigonometric model? If the model has nonlinear parameters in it, then you will need to use a nonlinear optimization. And depending on the model, to force it through a given point, that may require a nonlinear tool that can handle an equality constraint, so possibly fmincon. But sometimes, again, depending on the model, there may be simpler solutions. It all depends on the model.

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!