Help defining normalized SSE for problem based objective function

Hi all, I'm using problem-based solving to fit some coefficients in a differential equation to some experimental data that I'm reading in from an Excel file. Everything is working well. But I have some issues in that I'm working with multiple data sets simulaneously that span a wide range of responses and the sets at the high end of the distribution are forcing the fit at the expense of the datasets at the low end of the distribution just due to the fact that I'm using a sum-square-error objective function:

obj = sum(sum((myfcn - experimentalData{5:height(experimentalData),3:width(experimentalData)}).^2));

There are individual data sets in each column of the experimentalData table with different ranges. I want to scale the objective function calculation by the mean of these datasets.

I have calculated the mean for each one:

cAverage(i) = mean(experimentalData{5:height(experimentalData),i+2}); % there are two dummy columns before the data

What I want is a simple expression to divide the delta between the result of the calculated function and the experimental data by the average. So a matrix divided by a vector, by column.

When I type manually in the script prompt the right equation for this:

(myfcn - experimentalData{5:height(experimentalData),3:width(experimentalData)})./cAverage,

it appears to work, but when I edit the formula to

add this in and run the script, I'm getting an error for the argument dimensions must agree.

Can someone give me a clue how to do a normalized SSE objective function?

Thanks,

Tony

Answer 1

[Matt J]

add this in and run the script, I'm getting an error for the argument dimensions must agree.

Well, we have no idea what you've done to ensure that myfcn is the same size as your experimentalData block, but assuming you have, then,

 selection = experimentalData{5:height(experimentalData),3:width(experimentalData)};
 
 Error= ( myfcn-selection./mean(selection,1) );
 
 SSE=sum(  Error.^2,'all');

[Anthony Ozzello]

Thanks Matt.

the arrays are the same size, I have confirmed that. the normal SSE function works fine.

In the code you've given, we're dividing the experimental data by the mean of the experimental data if I read this right, but the function call returns un-normallized data. So I need to divide the difference between the function return and the experimental data by the mean. At least that's what I think I need to do. In this example, both myfcn and selection may be something like 13 x 3 with 13 data points from 3 different data series and they're both on the same scale and I believe I need to divide the difference between them by the mean of the experimental data series. Am I confused?

Thanks,

Tony

[Matt J]

That doesn't seem very different from what I proposed,

 Error = ( myfcn-selection)./mean(selection,1);

[Anthony Ozzello]

Thanks Matt,

Now, I'm relatively new to Matlab and I'm cobbling together pieces of code that I'm finding to make what I need work. For the problem-based setup. I have a single statement that defines the objective in terms of the function return and the experimental data and the function return is a call to my differential equation definition. The three calls look like this for the regular SSE calculation:

myfcn = fcn2optimexpr(@RtoODE,r,z,n,nComponents,tspan,y0, experimentalData, 'OutputSize',[height(experimentalData)-4,nComponents]);
obj = sum(sum((myfcn - experimentalData{5:height(experimentalData),3:width(experimentalData)}).^2));
prob = optimproblem("Objective",obj);

When I got with your solution to calculate the objective function, I won't be able to put the objective function definition on a single line, so I'll need to do it as a function call. I'm a little vague on how to set up the function call to the new objection function and reference myfcn in this case, would it look like this?

obj=myobj(@myfcn,r,z,n,nComponents,tspan,y0,experimentalData); ?

With the function 'myobj' then calculate and return the new SSE as you described? I can handle the writing of the function call. And then the call to myfcn is within the myobj function? Or do I somehow pass the output of myfcn to the myobj function call? I'm not sure.

Thanks and sorry for a potentially obvious question,

Tony

[Anthony Ozzello]

Hi Matt,

Just to check, I ran this code:

obj2 = myobj(myfcn,experimentalData);

function SSE = myobj(myfcn,experimentalData)
selection = experimentalData{5:height(experimentalData),3:width(experimentalData)};
Error = ( myfcn-selection)./mean(selection,1);
SSE=sum(  Error.^2,'all');
end

And I got the same error for the argument dimensions not matching. (myfcn-selection) is 13 x 2 and mean (selection,1) is 1 x 2.

Thanks,

Tony

[Matt J]

I can't tell what error you're seeing (because you didn't post it). I also can't advise how to set up the objective function because I can't see from your code which variables are the optimization variables and which are the fixed parameters.

I would recommend abandoning the problem-based framework for this, though. There's no real advantage to it when you don't have linear constraints. For this, I woudl just use lsqcurvefit().

Also, because this is an ODE optimization, you will probably need to be familiar with,

https://www.mathworks.com/help/optim/ug/optimizing-a-simulation-or-ordinary-differential-equation.html

[Anthony Ozzello]

Hi Matt,

Thanks for the suggestion. I have been through the fora on solving this. I have a system of differential equations and the basic outline for this on a different thread for using the method of lines to reduce a system of partial differential equations down to where I can use ode15s to address the stiff problem. I didn't want to overwhelm anyone with a ton of code. this is beyond what lsqcurvefit() can handle. Here is the full error message:

=====

Error using optim.internal.problemdef.ElementwiseOperator/checkIsValid

Argument dimensions 13-by-2 and 1-by-2 must agree.

Error in optim.internal.problemdef.Rdivide/checkIsValid

Error in optim.internal.problemdef.Rdivide.getRdivideOperator

Error in ./

Error in multicomponent_fitting_excel>myobj (line 199)

Error = ( myfcn-selection)./mean(selection,1);

Error in multicomponent_fitting_excel (line 73)

obj2 = myobj(myfcn,experimentalData); %trying to define normalized SSE

=====

I can attach the entire script if you like, but I didn't want to overload anyone. The whole script works great if I just use this for my objective function, as I stated above:

obj = sum(sum((myfcn - experimentalData{5:height(experimentalData),3:width(experimentalData)}).^2));

Thanks for being patient.

Tony

[Matt J]

this is beyond what lsqcurvefit() can handle

That's doubtful. SSE minimization is exactly what lsqcurvefit is for. In any case, there's no reaon to think the problem-based framework is going to do any better. It's just a wrapper for lsqcurvefit(0 and other solver-based methods.

In your code, you have this,

fcn2optimexpr(@RtoODE,r,z,n,nComponents,tspan,y0, experimentalData)

Which of these are variables that you are trying to optimize? One thing that's definitely wrong is that the list of variables is only supposed to include optimvars. Therefore, experimentalData should not be there.

[Anthony Ozzello]

Hi Matt,

r is the optimization vector consisting of several values. experimentalData is there because it gets passed to RtoODE which then calls the function that calculated my derivatives, so it's a passthrough. I didn't know any other way to get it down through subsequent calls to the final function that calculates the derivatives. Like this, I'm excluding some intervening code:

===

myfcn = fcn2optimexpr(@RtoODE,r,z,n,nComponents,tspan,y0, experimentalData);

obj = sum(sum((myfcn - experimentalData{5:height(experimentalData),3:width(experimentalData)}).^2));

prob = optimproblem("Objective",obj);

[rsol, sumsq] = solve(prob,r0,"Options",options);

function solpts = RtoODE(r,z,n,nComponents,tspan,y0, experimentalData)

[T, Y] = ode15s(@(t,y) Myfun(t,y,z,n,nComponents,r,experimentalData),tspan,y0,opts);

function DyDt=Myfun(t,y,z,n,nComponents,r,experimentalData)

*calculate the derivatives and return using data from, y,z,n,cComponents, and experimentalData using the method of lines in z and t

===

I hope that all makes sense.

I found this definition in the System Identification Toolbox which is exactly what I want to implement, just not sure how to code it:

'NMSE' fit(i)=||xref(:,i)−x(:,i)||^2/||xref(:,i)−mean(xref(:,i)||^2, fit is a row vector. then I just need to sum up all the values for fit(i).

[Anthony Ozzello]

the || indicates the 2-norm of the vector.

[Anthony Ozzello]

Hi Matt,

I got it working with the NMSE definition above for the objective function on a single line and everything is working well now. I appreciate the time you spend trying to help me figure this out.

Tony