Polynomial surface fitting problem with constraint y <=1

I'm trying to fit an emperical relationship between dependent variable y and two independent variables x1 & x2.

I tried to formulate it as an equation of form y = a*x1^b + c * x2^d where a, b, c, d are constants to be found. This is simple to do with lsqcurvefit. However, I also need to add a constraint that y <= 1, which is a property of the system. I tried looking into the Optimization Toolbox but can't seem to find a suitable solver.

Would really appreciate if someone could suggest either a way to solve this equation or a better way to formulate the relationship.

Below is the data points I have. I don't really know much about the relationship between y (capacity) and the two independent variables other than it decreases as the cumulative Ah discharged increases and is the highest at 25 degC.

Answer 1

[John D'Errico]

You have the model:

y = a*x1^b + c*x2^d

I'll assume that x1 is temperature, and x2 is current discharged.

You have told us it must be monotone decreasing with increasing discharge (x2). That just tells us that c*d must be negative, since the first derivative of that expression with respect to x2 is just

dy/dx2 = c*d*x2^(d-1)

With x2 positive, if that expression is negative, then c*d < 0 must apply.

But then you tell us the expression attains its maximum at T = x1 = 25. 25 is right in the middle of your test data. Can that expression have a maximum at some middle point in the domain? Again, form a partial derivative, this time with respect to x1.

dy/dx1 = a*b*x1^(b-1)

If it is maximal there, then you would have that expression == 0. And that means either x1=T=0, a=0, or b=0. Which is it? Any of those cases tells me your model makes no sense in terms of what you have told us.

I can stop right here, since your other constraint that the maximum must be 1 is irrelevant, IF the model makes no sense in context of what you have already told us. It would be easy to apply that constraint for this model, but still, why bother when the model is clearly not descriptive of what you have? Anyway, if the function must be maximal at T=x2=0, then all you need to do is constrain the value at the point T=0, C=25. But again, it is not relevant, since your function cannot encapsulate the curve shape you think must exist.

[Wenlin]

You are right, I chose the model only because it seemed like a simple curve to fit, but wasn't considering its actual shape. I guess instead of a*x1^b, I can use a*(x1-b)^c which would give me a peak along temperature = b?

[John D'Errico]

Um, within limits. That would be closer to being valid, assuming b was 25, so it would have a chance to have a maximum value there. But there is a serious problem even at that.

x1 varies from 15-35. If b was 25. then x1 - b will vary from -10 to 10. Correct? This term will be of the form a*(x1-b)^c.

where c would be presume to be some general real number. So try this:

(-10)^1.25

ans = -12.5743 - 12.5743i

What happens? Do you appreciate that a negative number, when raised to a non-integer power, will result in a complex number? As such you CANNOT write the expression as you want. Instead, consider this shape:

a*abs(x1 - b)^c

Now, if c is a positive real number, and b is 25, then it will be closer to what you want. a must be negative.

fplot(@(x) -1*abs(x-25).^1.25,[15 35])

Depending on the power, this curve will fall off at varying rates.

So you might hop to try a model like this:

y = a*abs(x1 - 25).^b + c*x2^d

Does that work? NO. Still not yet.

a must be negative there. b must be positive. And we have the additional constraint that c*d must be negative. We know that much so far.

Now think about what value this form takes on when T=x1= 25, and C=x2=0? That MUST be the maximum value. And it should be something on the order of 1, right? So try it!

y(25,0) = a*abs(25 - 25).^b + c*0^d = 0

So the maximum possible value that relationship can take on will be ZERO. You are still missing one term, i.e., a constant. That constant will allow the entire surface to translate (or shift if you prefer) up or down along the z axis. The form you will need must be more like this:

y = a*abs(x1 - 25).^b + c*x2.^d + e

there are now 5 unknown parameters in the model. e is the maximum value, the value this relationship will take on at T=25, and C=0. So you now have one more constraint, a bound constraint that e<=1.

With those changes, your function will now have the desired shape. For example, with some random set of parameters, it might look like this:

fsurf(@(x1,x2) -.001*abs(x1-25).^2 - .00002*x2.^1.25 + 1,[15 35 0 1200])
xlabel 'Temperature'
ylabel 'Ah Discharged'

As you can see, it is a decreasing function of discharge, with a peak at T=25 and C=0. The maximum value is purely controlled by e, so as long as e is no larger than 1, you are fine. It falls off in both directions as T moves away form 25.

[Wenlin]

This makes a lot of sense and I was able to fit my data perfectly. Thank you so much!

[John D'Errico]

Let me add, that this is a common way one builds such an empirical model. You start with a fundamental shape that may be appropriate for the model in one variable. Then add to it extra terms to add in the necessary behavior. Eventually you find a reasonable model for the process, one with terms that can then be estimated from the data.

[Wenlin]

Yeah that makes sense, thank you for guiding me through it!

Answer 2

[Torsten]

If N is the number of (y12,x1,x2) data, use fmincon instead of lsqcurvefit and set the constraints

a*x1_i^b + c*x2_i^d - 1 <= 0 (i=1,2,...,N)

in the function "nonlcon".