Optimization single nonlinear equation to fit predicted data to actual data

I am trying to come up with an equation that will use parameters of a given object to predict how it will function. I have actual performance data (call it ydata) and 3 parameters of 11 different objects as my data set. I start with a guess (call it xdata) for each object and am trying to optimize an equation of the format a * p1^w * p2^y * p3^z = xdata such that xdata is as close as possible to y data for each point (in this case, p1, p2, p3 are the known parameters and w, y, z are the variables that need to be changed to get xdata to match ydata and a is just a constant). What I want to end up with is one value for w, one value for y, and one value for z that optimizes all 11 data points at once, essentially making sure the R^2 value comparing xdata and ydata of all 11 points is >=0.9 but I cannot figure out a way to do this successfully. It seems like no matter what method I try (lsqcurvefit, surrogateopt, solve, fmincon, etc) runs into some sort of error or will not work with the problem I am trying to solve. For most of the optimization methods that look to minimize, I have tried setting up the equation like so: a * p1^w * p2^y * p3^z - xdata = 0 and setting the constraint that the function value must be greater than 0 but I can't figure out a way to constrain the function value itself, it seems only the variables you are trying to change (w,y,z) can be constrained.

Any thoughts, even a general direction, would be super helpful at this point.

[othman warde]

[١٥/٢, ١٢:٣٢ ص] .othman Alibory: % Define the function to be optimized fun = @(p, xdata) p(1) * xdata(:, 1).^p(2) .* xdata(:, 2).^p(3) .* xdata(:, 3).^p(4);

% Define the initial guess for the parameters p0 = [1, 1, 1, 1];

% Load the data load('data.mat'); % replace with the name of your data file xdata = data(:, 1:3); ydata = data(:, 4);

% Define the nonlinear least squares problem problem = struct; problem.objective = @(p) ydata - fun(p, xdata); problem.x0 = p0; problem.lb = [0, -Inf, -Inf, -Inf]; % Lower bounds on the parameters problem.ub = [Inf, Inf, Inf, Inf]; % Upper bounds on the parameters

% Solve the nonlinear least squares problem options = optimoptions('lsqnonlin', 'Algorithm', 'trust-region-reflective', 'Display', 'iter'); [p, resnorm, residual, exitflag, output, lambda, jacobian] = lsqnonlin(problem);

% Extract the optimized parameters a = p(1); w = p(2); y = p(3); z = p(4);

% Evaluate the model with the optimized parameters xdata_fit = xdata; ydata_fit = fun(p, xdata);

% Compute the R^2 value R2 = 1 - sum((ydata - ydata_fit).^2) / sum((ydata - mean(ydata)).^2);

% Plot the results figure; scatter(ydata, ydata_fit); xlabel('Actual Performance Data'); ylabel('Predicted Performance Data'); title(sprintf('R^2 = %.2f', R2)); [١٥/٢, ١٢:٣٢ ص] .othman Alibory: This code defines a function fun that computes the predicted performance data given the parameters and the input data, and sets up a nonlinear least squares problem using the lsqnonlin function. The initial guess for the parameters is set to p0 = [1, 1, 1, 1], but you can modify this as needed. The lower and upper bounds on the parameters are set to [0, -Inf, -Inf, -Inf] and [Inf, Inf, Inf, Inf], respectively, to ensure that the parameters are positive. The lsqnonlin function then searches for the values of the parameters that minimize the sum of squared residuals between the observed and predicted values.

After finding the optimized parameters, the code evaluates the model with these parameters and computes the R^2 value to evaluate the goodness of fit. Finally, the code plots the actual versus predicted performance data to visualize the results.

Note that this is just an example, and you may need to modify the code to fit your specific problem. In particular, you should replace data.mat with the name of your data file and modify the input and output variables of the fun function to match your parameter and data format.

Answer 1

[Torsten]

"nonlcon" as the nonlinear constraint function from "fmincon" can constrain a*p1(i)^w * p2(i)^y * p3(i)^z as well, if this is what you want.

But maybe - to get a good initial guess for a,w,y and z - you should first take log of your model and thus fit a linear model

log(ydata) = log(a) + w*log(p1) + y*log(p2) + z*log(p3)

where ydata are your (measured) performance data and p1, p2 and p3 are the factors that influence this performance.

[Mckenzie]

Taking the log of the equation is an exellent suggestion, I will give that a go. As for using fmincon, what would that set up look like? I have tried the following:

fun = @(x) 5.4 * 5^x(1) * 18^x(2) * 0.6^x(3);
A = [1,1,1];
b = [0];
Aeq = [];
beq = [];
lb = [0.2,0.2,0.2];
ub = [1,1,1];
x0 = [0.5,0.5,0.5];
nonlcon = @ineq;
[x,fval] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)

where @ineq is a saved function that is:

function [c,ceq] = ineq(x)
c = 5.4 * 5^x(1) * 18^x(2) * 0.6^x(3) - 87.1 - 87.1*0.05;
ceq = [];

Note here that a = 5.4, p1 = 5, p2 = 18 and p3 = 0.6, also the actual value (ydata) that I am trying to guess is 87.1 and I set up the inequality to be a * p1^w * p2^y * p3^z <= actual+5%. This doesn't work the way I want, the value for fval is usually -70 where it should be something close to 87.1

Am I doing something incorrectly? I am not sure where to put the constraint to say fval (also x data) needs to be >0, also i have not attempted to iterate in this edition, I just want to make sure it's doing what I intend.

[Torsten]

You set the constraint that x(1) + x(2) + x(3) <= 0, but you also set x(1),x(2),x(3) >= 0.2. This cannot work.

[Mckenzie]

Ok, I have deleted those conditions and are getting answers that make sense, thanks! Now for iteration with function handles. I am trying the following:

for n = 1:size(num,1)
    fun = @(x) coeff(n) * num(n,3)^x(1) * num(n,4)^x(2) * num(n,5)^x(3);
    A = [];
    b = [];
    Aeq = [];
    beq = [];
    lb = [0.5,0.2,0.5];
    ub = [0.9,0.6,0.9];
    x0 = [0.6,0.6,0.6];
    nonlcon = @ineq;
    [x,fval] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)
end

where a, p1,p2,p3 are pulling from matricies coeff and num. I am getting an error with the nonlcon line. I have updated my ineq function to be:

function [c,ceq] = ineq(x)
c = [];
ceq = coeff(n) * num(n,3)^x(1) * num(n,4)^x(2) * num(n,5)^x(3) - num(n,1);

(changing from the inequality method to the equality one because I got better results). It doesn't like that I have the "n" in the function handle but I need that function to update with each interation. How can I achieve this?

[Torsten]

As far as I understood your question, you want to determine x(1),x(2),x(3) and x(4) such that

sum_{i=1}^{i=size(num,1)} (x(1) * num(i,3)^x(2) * num(i,4)^x(3) * num(i,5)^x(4) - num(i,1))^2

is minimized (maybe with some constraints on the parameters and the function values).

This doesn't correspond to your above setup. But maybe I got something wrong.

[Mckenzie]

I want to find x(1),x(2),x(3) such that the equation coeff(n) * num(n,3)^x(1) * num(n,4)^x(2) * num(n,5)^x(3) - num(n,1) is minimize but positive (so as close to zero as possible). I can do this for one value of my data set by hardcoding the matrix location of the value I working on, but where I run into problems is iteration through the entire data set. Here is my most recent attempt:

for n = 1:size(num,1)
    fun = @(x) coeff(n) * num(n,3)^x(1) * num(n,4)^x(2) * num(n,5)^x(3) - num(n,1);
    A = [];
    b = [];
    Aeq = [];
    beq = [];
    lb = [0.2,0.2,0.2];
    ub = [0.9,0.9,0.9];
    x0 = [0.6,0.6,0.6];
    nonlcon = @(x)ineq(x(n));
    [x,fval] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)
    
end

function [c,ceq] = ineq(x)
[num,txt,raw] = xlsread('test.xlsx', 'Data');
coeff = num(:,2);
c = coeff(n) * num(n,3).^x(1) * num(n,4).^x(2) * num(n,5).^x(3) - num(n,1) - num(n,1)*0.5;
ceq = coeff(n) * num(n,3)^x(1) * num(n,4)^x(2) * num(n,5)^x(3) - num(n,1);
end

But I get an error saying "undefine function or variable 'n' within the funtion ineq" am I approching this wrong?

[Torsten]

As far as I understand your problem, it has to be posed this way:

[num,txt,raw] = xlsread('test.xlsx', 'Data');
fun = @(p) sum((num(:,2).*num(:,3).^p(1) .* num(:,4).^p(2) .* num(:,5).^p(3) - num(:,1)).^2);
A = [];
b = [];
Aeq = [];
beq = [];
lb = [0.2,0.2,0.2];
ub = [0.9,0.9,0.9];
p0 = [0.6,0.6,0.6];
[p,fval] = fmincon(fun,p0,A,b,Aeq,beq,lb,ub,@(p)ineq(p,num))

function [c,ceq] = ineq(p)
  c = -(num(:,2) .* num(:,3).^p(1) .* num(:,4).^p(2) .* num(:,5).^p(3) - num(:,1)) ;
  ceq = [];
end

[Mckenzie]

The goal is to have it come up with 11 different solutions (one for each data point) for x1, x2, and x3 which is why I was looking into the for loop, so a simple summation would not work in this context.

[Torsten]

You want different parameters for each measurement ? But then you have 11x1 equation for 11x3 unknowns which will give you -eh, bullshit.

The parameters will be almost arbitrary because you can always fix 2 parameters and find a remaining one that makes

num(:,2) .* num(:,3).^p(1) .* num(:,4).^p(2) .* num(:,5).^p(3) - num(:,1) == 0

Take e.g.

p(1) = p(2) = 0 

and

 p(3) = log(num(:,1)/num(:,2)) / log(num(:,5))