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

1 次查看（过去 30 天）

显示更早的评论

Dursman Mchabe 2019-1-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/439078-unable-to-perform-assignment-because-the-left-and-right-sides-have-a-different-number-of-elements

评论： Dursman Mchabe 2019-1-11

DynamicOptimization10Jan2019.m

在 MATLAB Online 中打开

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

16 个评论
显示 14更早的评论隐藏 14更早的评论

Dursman Mchabe 2019-1-10

在 MATLAB Online 中打开

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?

Dursman Mchabe 2019-1-10

在 MATLAB Online 中打开

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.

Dursman Mchabe 2019-1-10

The initial conditions are:

x0 = 9.54992586;

y0 = 19.89;

z0 = 0;

w0 = 0;

Dursman Mchabe 2019-1-11

在 MATLAB Online 中打开

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