Unable to perform assignment because the left and right sides have a different number of elements.

Hi all

When I run the attached code, I get the error message:

Unable to perform assignment because the left and right sides have a different number of elements.

Error in DynamicOptimization10Jan2019>myModel (line 90)
dxdt(1)= -2.85.*(1/((1/0.00639039)+(1/k)).*12.54.*x.*x.*(1-(2.69636.*x)/(1+(2.69636.*x))));

Error in DynamicOptimization10Jan2019>@(t,x)myModel(t,x,k) (line 81)
[t,x] = ode15s(@(t,x)myModel(t,x,k),time,x0);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode15s (line 150)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in DynamicOptimization10Jan2019>myObj (line 81)
[t,x] = ode15s(@(t,x)myModel(t,x,k),time,x0);

Error in fminsearch (line 200)
fv(:,1) = funfcn(x,varargin{:});

Error in DynamicOptimization10Jan2019 (line 3)
k = fminsearch(@myObj,k0);
 

I know the error I am getting implies that I am trying to put more elements into less elements, or vice versa. However, I do not know how to make elements equal.

please help.

Kind Regards

[KSSV]

YOur output in the line:

dxdt(1)= -2.85.*(1/((1/0.00639039)+(1/k)).*12.54.*x.*x.*(1-(2.69636.*x)/(1+(2.69636.*x))));

is a matrix..and you are trying to save it as array. You need to rethink on your code.

[Dursman Mchabe]

Thanks a lot such a quick response. I will rethink. My first move will be to replace .* with *.

[Dursman Mchabe]

Oops, it gives another error:

Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform
elementwise multiplication, use '.*'.

Error in DynamicOptimization10Jan2019>myModel (line 89)
dxdt(1)= -2.85*(1/((1/0.00639039)+(1/k))*12.54*x*x*(1-(2.69636*x)/(1+(2.69636*x))));

I must rethink.

[Walter Roberson]

In particular, inside myModel, x is going to be the same size as x0 .

[Dursman Mchabe]

Thanks a lot for the comment Walter. The challenge is that I have 4 ODEs dxdt(1) ... dxdt(4), with only 2 varibles, x(1) and x(2). Perhaps I should try to pass x as x(1) and x(2).

[Dursman Mchabe]

Oops, it gives:

Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 1-by-4.

Error in fminsearch (line 200)
fv(:,1) = funfcn(x,varargin{:});

Error in DynamicOptimization10Jan2019 (line 3)
k = fminsearch(@myObj,k0);

How can I change the left to a 1-by-4?

[Walter Roberson]

You cannot change the left to 1 x 4. You are using myObj as the objective function for fminsearch, and objective functions for fminsearch must always return scalars. Returning 1 x 4 there would be like trying to simultaneously minimize 4 things, but without instructions about the relative values of decreases between the four parts (because you could easily encounter a situation where a small increase in one of the four is needed to permit a large decrease in one or more of the other three.)

If you do want simultaneous minimization of four outputs, then you should look at gamultiobj and examine the pareto front.

[Dursman Mchabe]

I want to try and keep things simple. When I hypothetically solve just one ODE, the code works fine. See below:

%% solve with fminsearch
k0 = 59563518.8216;
k = fminsearch(@myObj,k0);
disp(['fminsearch: k = ' num2str(k)]);
mySim(k);

%% solve with fmincon
LB = 20;
UB = 100;
k = fmincon(@myObj,k0,[],[],[],[],LB,UB); 
disp(['fmincon: k = ' num2str(k)]);
mySim(k);

%% solve graphically
n = 100;
k = linspace(4e8,5e9,n);
for i = 1:n,
    obj(i) = myObj(k(i));
end
figure(2)
plot(k,obj,'b--','LineWidth',3);
xlabel('k value')
ylabel('Objective Value')

%% confidence interval
%(S(k) - S(K)) / S(K) <= p / (n-p) * F(p,n-p,1-alpha)
p = 1; % number of parameters
np = 30; % number of data points
alpha = 0.05; % alpha, confidence interval
rhs = p / (np-p) * finv(1-alpha,p,(np-p));
ub = min(obj) * rhs + min(obj);
hold on
plot([min(k) max(k)],[ub ub],'r-','LineWidth',3);
legend('Objective','Confidence Limit')
% axis([0.06 0.09 0 1])
% axis([0.07 0.09 0 1])

axis([3e8 5e8 0 7e9])
axis([4e8 5e8 0 7e9])

function obj = myObj(k)

% measured values

time = [1
3
5
7
9
11
13
15
17
19
21
23
25
27
29
31
];


y = [9.54992586
5.623413252
2.884031503
0.489778819
0.048977882
0.046773514
0.032359366
0.022908677
0.016595869
0.014125375
0.011748976
0.01 
0.00851138
0.007585776
0.00691831
0.006309573
];

% predicted values

x0 = 9.54992586;
[t,x] = ode15s(@(t,x)myModel(t,x,k),time,x0);

% compute sum of squared errors
obj = sum((x-y).^2);
end

function dxdt = myModel(t,x,k)

dxdt= -2.85*(1/((1/0.00639039)+(1/k))*12.54*x*x*(1-(2.69636*x)/(1+(2.69636*x))));
end


function [] = mySim(k0)

% Predicted values
x0 = 9.54992586;
[t,x] = ode15s(@(t,x)myModel(t,x,k0),[0 31],x0);

% Actual values
time = [1
3
5
7
9
11
13
15
17
19
21
23
25
27
29
31
];


y = [9.54992586
5.623413252
2.884031503
0.489778819
0.048977882
0.046773514
0.032359366
0.022908677
0.016595869
0.014125375
0.011748976
0.01 
0.00851138
0.007585776
0.00691831
0.006309573
];



figure(1)
hold off
plot(t,x)
hold on;
plot(time,y,'ro')
xlabel('Time')
ylabel('Value')
legend('Predicted','Actual')
end

Is there a better way of applying this code on more than one ODEs?

[Walter Roberson]

It is pretty late where I am, but at the moment I get the impression that this is effectively a boundary value problem, trying to find the k that leads to a particular outcome.

I am also suspecting that your model has a closed form solution if approached analytically, so I not certain you need the ode15s call. However, I am too tired to figure out what your original equations must have been.

[Dursman Mchabe]

I'm sorry to keep you up. I truly appreciate your help. Where I am it is 11:32 am.

You are absolutely right. I am looking for a k that will make model calculations equal measured values.

[Dursman Mchabe]

When I take lot of gueses manually, I find k0 = 46.83370 to work see below:

function [] = mySim(k0)
k0 = 46.83370;
% Predicted values
x0 = [9.54992586;19.89;0;0];
[t,x] = ode15s(@(t,x)myModel(t,x,k0),[0 31],x0);

% Actual values
time = [1
3
5
7
9
11
13
15
17
19
21
23
25
27
29
31
];


y = [9.54992586	19.89	0	0
5.623413252	18.51227628	1.377723722	0.000176583
2.884031503	17.5510897	2.338910301	0.000584451
0.489778819	16.71100104	3.178998962	0.004671909
0.048977882	16.55633404	3.333665957	0.048352564
0.046773514	16.55556058	3.33443942	0.050608503
0.032359366	16.55050298	3.339497016	0.072770512
0.022908677	16.54718695	3.342813047	0.101976407
0.016595869	16.54497193	3.345028067	0.139244055
0.014125375	16.54410509	3.345894907	0.16245708
0.011748976	16.54327127	3.346728731	0.193464978
0.01        16.54265759	3.347342407	0.225067572
0.00851138	16.54213527	3.34786473	0.26139843
0.007585776	16.5418105	3.348189503	0.290553969
0.00691831	16.5415763	3.348423702	0.315962812
0.006309573	16.54136271	3.348637294	0.34334241];


figure(1)
hold off
plot(t,x)
hold on;
plot(time,y,'v')
xlim([0 35])
ylim([-1 20])
xlabel('Time (sec)')
ylabel('Concentration (mol/m^3)')
legend('Exp-H^+', 'Exp-CaCO_3', 'Exp-Ca^2+', 'Exp-HCO^-_3','Mod-H^+', 'Mod-CaCO_3', 'Mod-Ca^2+', 'Mod-HCO^-_3', 'Location','best')
end

function dxdt = myModel(t,x,k)

dxdt=zeros(4,1);
dxdt(1)= -2.85*(1/((1/0.00639039)+(1/k))*12.54*x(2)*x(1)*(1-(2.69636*x(1))/(1+(2.69636*x(1)))));
dxdt(2)= -(1/((1/0.00639039)+(1/k))*12.54*x(2)*x(1)*(1-(2.69636*x(1))/(1+(2.69636*x(1)))));
dxdt(3)= (1/((1/0.00639039)+(1/k))*12.54*x(2)*x(1)*(1-(2.69636*x(1))/(1+(2.69636*x(1)))));
dxdt(4)= (1/((1/0.00639039)+(1/k))*12.54*x(2)*x(1)*(1-(2.69636*x(1))/(1+(2.69636*x(1)))));

end

How ever I need to use an objective function for a lot other similar codes.

[Walter Roberson]

can you post the equation form of the differentiation equationss?

You should also look at bvp4c and related functions .

[Dursman Mchabe]

[Dursman Mchabe]

Sorry for the delay.

[Dursman Mchabe]

The initial conditions are:

x0 = 9.54992586;

y0 = 19.89;

z0 = 0;

w0 = 0;

[Dursman Mchabe]

I have found partial success when using lsqcurvefit. It can estimate k. Yet I don't know how to do sensivity analysis and to determine confidence interval without the objective function. Please see the code below.

function KineticModelat30deg

    function C=KineticModel(k,t)
    c0=[9.54992586;19.89;0;0];
    [T,Cv]=ode45(@ODEs,t,c0);
%
    function dC=ODEs(t,c)
    dcdt=zeros(4,1);
    dcdt(1)= -2.85*(1/((1/0.00639039)+(1/k(1)))*12.54*c(2)*c(1)*(1-(k(2)*c(1))/(1+(k(2)*c(1)))));
    dcdt(2)= -(1/((1/0.00639039)+(1/k(1)))*12.54*c(2)*c(1)*(1-(k(2)*c(1))/(1+(k(2)*c(1)))));
    dcdt(3)= (1/((1/0.00639039)+(1/k(1)))*12.54*c(2)*c(1)*(1-(k(2)*c(1))/(1+(k(2)*c(1)))));
    dcdt(4)= (1/((1/0.00639039)+(1/k(1)))*12.54*c(2)*c(1)*(1-(k(2)*c(1))/(1+(k(2)*c(1)))));
    dC=dcdt;
    end
C=Cv;
end

t=[1
3
5
7
9
11
13
15
17
19
21
23
25
27
29
31
];

c=[9.54992586	19.89	0	0
5.623413252	18.51227628	1.377723722	0.000176583
2.884031503	17.5510897	2.338910301	0.000584451
0.489778819	16.71100104	3.178998962	0.004671909
0.048977882	16.55633404	3.333665957	0.048352564
0.046773514	16.55556058	3.33443942	0.050608503
0.032359366	16.55050298	3.339497016	0.072770512
0.022908677	16.54718695	3.342813047	0.101976407
0.016595869	16.54497193	3.345028067	0.139244055
0.014125375	16.54410509	3.345894907	0.16245708
0.011748976	16.54327127	3.346728731	0.193464978
0.01	16.54265759	3.347342407	0.225067572
0.00851138	16.54213527	3.34786473	0.26139843
0.007585776	16.5418105	3.348189503	0.290553969
0.00691831	16.5415763	3.348423702	0.315962812
0.006309573	16.54136271	3.348637294	0.34334241
];

k0=[0.1;50];

[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@KineticModel,k0,t,c);

fprintf(1,'\tRate and Adsorption constants:\n')
for i = 1:length(k)
    fprintf(1, '\t\tk(%d) = %8.5f\n', i, k(i))
end

tv = linspace(min(t), max(t));
Cfit = KineticModel(k, tv);

figure(1)
plot(t, c, 'v')
hold on
plot(tv, Cfit);
hold off
xlabel('Time (sec)')
ylabel('Concentration (mol/m^3)')
legend('Exp-H^+', 'Exp-CaCO_3', 'Exp-Ca^2+', 'Exp-HCO^-_3','Mod-H^+', 'Mod-CaCO_3', 'Mod-Ca^2+', 'Mod-HCO^-_3', 'Location','E')

end