ODE45 and lsqcurvefit - Function value and YDATA sizes are not equal.

Dear MATLAB Community,

I have an issue with parameter optimization using ODE45 and lsqcurvefit. I've noticed that many users had simillar problem but each case has to be dealt separatly.

I would like to find the values of two parameters in system of ODEs. The equations are a little bit complicated (it is a QMOM opproach to Population Balance Equations). If you would like to see the whole code please let me know.

The script looks as follows:

function OPT

function M=PBE(parametry,t)
 
    m0(1) = 1;
    m0(2) = 1;
    m0(3) = 1;
    m0(4) = 1;
    m0(5) = 1;
    m0(6) = 1;
    
    N = 3;
     
    opts = odeset('RelTol',1e-3,'AbsTol',1e-6,'MaxStep',0.001);
    [T,Mv]=ode45(@(t,m)dmdt(t,m,N,parametry),t,m0,opts);
        
    M=Mv;
    
end

t=[1
    3]

m=[0.0559 0.0937 0.2475  1.0000 5.0708 29.0395
   0.0058 0.0264 0.1562 1.0000 6.5717 43.8660]

parametry0=[1.00;0.1]

PBE(parametry0,[0 1])


[parametry]=lsqcurvefit(@PBE,parametry0,t,m)

where dmdt is a set of ODEs describing of population moments change in time. dmdt is the function of time, moments, numbers of nodes N (simply speaking the ODE system consist of 2N equations), and two constant parameters parametry(1) and parametry(2).

I am trying to fit the functions to two sets of points for 1st and 3rd second (in the script correct values are present for the guess parameters value parametry0(1)=1 and parametry0(2)=0.1 just to check if the code is working). I am receiving the following result:

t =

     1
     3


m =

    0.0559    0.0937    0.2475    1.0000    5.0708   29.0395
    0.0058    0.0264    0.1562    1.0000    6.5717   43.8660


parametry0 =

    1.0000
    0.1000


ans =

    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000
    0.9990    0.9993    0.9996    1.0000    1.0005    1.0012
    0.9980    0.9985    0.9992    1.0000    1.0010    1.0023
    
                            ...(multiple results lines)
    
    0.0560    0.0937    0.2476    1.0000    5.0695   29.0305
    0.0559    0.0937    0.2476    1.0000    5.0702   29.0350
    0.0559    0.0937    0.2475    1.0000    5.0708   29.0395

Error using lsqcurvefit (line 271)
Function value and YDATA sizes are not equal.

Error in OPT (line 36)
[parametry]=lsqcurvefit(@PBE,parametry0,t,m)

So, as you can see the PBE function works correctly. It generates 6 columns of data. The time matrix has 2 rows. The YDATA matrix has 2 rows and 6 columns. The initial guess matrix has 2 rows. So what am I doing wrong?

Moreover, I have one additional question. How does one inputs the initial time value in such script? According to manual tspan=[initial_time final_time]. Yet here, one only enters the final times. Is it assumed by default that t_initial=0? Is it possible to change this assumption?

Best regards,

Marcin

[Star Strider]

Post the relevant parts of the code if you want help with this problem.

For a version of that problem that works correctly, see: Parameter Estimation for a System of Differential Equations

Since the error is:

Error using lsqcurvefit (line 271)
Function value and YDATA sizes are not equal.

and the correct lsqcurvefit syntax is:

x = lsqcurvefit(fun,x0,xdata,ydata)

it is obvious that whatever ‘m0’ is (that we do not have) does not have the same row and column dimensions of the output of ‘PBE’.

You need to be certain that it does, in order to do the regression you want.

Answer 1

[Marcin Lemanowicz]

Thank you for your answer. I am attaching the model files.

Some explenation:

dmdt.m represents the system of ODEs. It may be written as

I am sure that it could be written more professional but I am a beginner in the matlab environment. Since the first row and column in matlab have index 1 the 0 order moment is written as the 1st equation etc. The first term on the right side of the equation represents the birth of particles due aggregation (B_A.m file). Beta(i,j) represents the aggregation kernel (Beta_A.m), whereas w and L represents weights and abscissas received using Product Difference algoritm (PD.m). The second term represents the death due aggregation (D_A.m). The third term represents the birth due breakage (B_B.m). Here a(i) stands for the breakage kernel (Beta_B.m) and b(i) is the fragment distribution function (gama.m).

The mathematical model works perfectly when it is used directly. It calculates the evolution of the moments in time. However if I would like to find the parameters the issue described above appears.

Best regards,

Marcin

Answer 2

[Marcin Lemanowicz]

Dear MATLAB Community,

I've solved the problem. Simply speaking for the initial guessed parameters and the target parametrs the mathematical model gave real numbers. However, during the optimisation process a specific set of parametres was found which generated complex numbers. Therefore the size error occured. I had to decrease the tolerance value significantly for the ode45 function.

However I've encountered another issue. When I'm trying to find the solution using lsqcurvefit I am obtaining the following result:

It's quite good but not perfect.

Therefore I've tried to find the solution using fminsearch function. I've constructed a target function which is a sum of squares of differences between given and calculated values. The results are far more better:

Could you please help me with that issue? The initial guessed values are identical. How to enhance the solution in the first case?

I am sending all files. The file structure is the same as discussed earliel. CALCULATIONS presents simple calculations without optimisation procedure. OPT_LSQCURVEFIT presents the first case whereas th OPT_FMINSEARCH presents the second case.

[Alan Weiss]

I would ask/suggest several things:

• If plotting the logarithms of the results is the right thing, then perhaps you should optimize the logarithms of your responses.
• I suspect that you should put bounds on some or all of your variables, such as a lower bound of 0 on some control variables. This setting can make lsqcurvefit much more robust.
• Do you calculate the sum of squared differences in your objective function code, or do you pass the vector of response values? Despite the documentation warning not to calculate the sum of squares inside your code, I have found many times that is what people do, and it makes the results much poorer.
• Diid you try starting lsqcurvefit from several initial points?
• Maybe you need to set larger-than-default finite differences as suggested in Optimizing a Simulation or Ordinary Differential Equation (ODE)

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation