The error 'Matrix dimensions must agree' during optimization.How can it be solved?

I have the following code where 'Ymeas' is an array of measured values and must be fitted with estimated values-'Yest' meanwhile optimizing K1 using an ode.The objective function is the sum of squares of difference b/w Ymeas and Yest. Ymeas and Yest are 11X2 arrays. It optimizes for a while and then the following error is shown:

'Matrix dimensions must agree.'

 function optimize

 [K1] = deal(3.5); 

 Ymeas = [
         1  1

        0.8696  0.8696
        0.7692  0.7692
        0.6897  0.6897
        0.625  0.625
        0.5714  0.5714
        0.5263  0.5263
        0.4878  0.4878
        0.4545  0.4545
        0.4255  0.4255
        0.4  0.4     ];

 T=[0:10:100];

 A0=Ymeas(1,1);
 B0=Ymeas(1,2);

 figure;
 plot(T,Ymeas);
 hold on;
 h = plot(T,Ymeas,'k','linewidth',2);
 minERR = Inf;
 opts = optimset('fminunc');
 opts.LargeScale = 'off';
 Xest = fminunc(@(X)objfun(X),[0],opts);
 Xest = num2cell(Xest);
 [K1] = deal(Xest{:}),

           function ERR = objfun(X);
               X = num2cell(X);
              [K1] = deal(X{:});
   [T,Yest] = ode45(@(t,y)[-K1*y(1)*y(2); 
                           -K1*y(1)*y(2)],[T],[A0,B0]);
              ERR = sum((Ymeas(:) - Yest(:)).^2);
              if ERR < minERR
                  minERR = ERR;
                  for n = 1:2; set(h(n),'Ydata',Yest(:,n)); end
                  drawnow;
              end;
          end;
  end

[Walter Roberson]

What are you expecting from your line

[K1] = deal(X{:});

The only point of that I can see at the moment is to issue an error message if X was not a scalar.

[Nitin Samuel]

oh..but even without that , how can there be a problem with the matrix dimensions? Dimensions of both Ymeas and Yest are constant and only the values of Yest change...So why exactly is there an issue with the dimensions? Thanks...

Answer 1

[Teja Muppirala]

The problem occurs when K1 is negative. The ODE solver blows up before it can integrate to all your time values. You should specify a constraint that forces K1 to be >= 0. For example, you might use FMINCON instead of FMINUNC.

[Nitin Samuel]

Can you send me a sample code with only the 2-3 statements of code for fmincon I should use to replace this one? Im new to optimization and it would be really helpful.

Thank you Mr.Teja.

[Teja Muppirala]

Xest = fmincon(@(X)objfun(X),[0],[],[],[],[],0,[],[],opts);

See the help for FMINCON

[Nitin Samuel]

Thank you so much..it works perfectly...

[Nitin Samuel]

For K1 and K2, what should I vary in that statement you put?I tried adding another zero as in

Xest = fmincon(@(X)objfun(X),[0 0],[],[],[],[],0,[],[],opts);

But it doesnt work..It would be of great help.Thank you..

[Alan Weiss]

Take a look at the syntax for bounds. You should write

 Xest = fmincon(@(X)objfun(X),[0,0],[],[],[],[],[0,0],[],[],opts);

Or pick a better point than [0,0] for x0.

Alan Weiss

MATLAB mathematical toolbox documentation

[Nitin Samuel]

Thank you :)