Problems using ODE45 using fmincon command for least square optimisation problem
1 次查看(过去 30 天)
显示 更早的评论
Script file
function chisquare=myObjective(P)
expdata=load('ru.txt');
straindata=expdata(:,2);
Timedata=expdata(:,1);
a=0.5179;
b=1.1003E-3;
c=0.4576;
d=1;
sigma=180;
hstress=3463.46;
Hstar=0.7846;
Kc=0.1137;
tspan=0:0.01:40;
initial_conditions=[0;0;0;0];
function dy=pair1(t,y,a,b,d,hstress,sigma,Hstar,kc,c) %function definition
dy=zeros(4,1);
dy(1)=a*sinh((b*(sigma^d)*(1-y(2)))/((1-y(3))*(1-y(4)))) ;
dy(2)=((hstress*(a*sinh(b*(sigma^d)*(1-y(2)))/(1-y(3))*(1-y(4))))/(sigma^d))*(1-y(2)/Hstar);
dy(3)=(kc/3)*power((1-y(3)),4);
dy(4)=c*y(1);
end
end
[t,y]=ode45(@(t,y)pair(t,y,a,b,d,hstress,sigma,Hstar,Kc,c),tspan,initial_conditions) %function call
plot(t,y(:,1))
diffsquare(:,i)=((expdata(i,2)-y(1)).^2)/n;
end
chisquare=sum(diffsquare);
global c;
c=(chisquare);
Second File
P0=[7E-6; 5E-6; 0.5; 0.43];
lb=[1E-6; 5E-8; 0; 0.42];
ub=[3E-5; 9E-6; 1; 0.5];
P=fmincon(@myObjective,P0,[],[],[],[],lb,ub);
plot(t,Straindata,'+-','Linewidth',1.5,'color','red');
hold on
plot(t,y(1),'--','Linewidth',1.5,'color','blue');
xlabel('Time');
ylabel('Strain');
Ending with the error
Error using feval
Error: File: myObjective.m Line: 24 Column: 1
This statement is not inside any function.
(It follows the END that terminates the definition of the function "myObjective".)
Error in fmincon (line 534)
initVals.f = feval(funfcn{3},X,varargin{:});
Error in Untitled8 (line 4)
P=fmincon(@myObjective,P0,[],[],[],[],lb,ub);
Caused by:
Failure in initial user-supplied objective function evaluation. FMINCON cannot continue.
I suspect I know a bit of what's going wrong.. yet I'm at a loss for how to fix it.. Any help is much appreciated
0 个评论
回答(2 个)
Bjorn Gustavsson
2019-7-13
Your first function is messed up with code after the end that closes the function myObjective, try this change
function chisquare=myObjective(P,straindata,Timedata)
% expdata=load('ru.txt'); % Better to send in straindata
% straindata=expdata(:,2); % and Timedata as parameters than
% Timedata=expdata(:,1); % to load them every time you call the function
a = 0.5179;
b = 1.1003E-3;
c = 0.4576;
d = 1;
sigma = 180;
hstress = 3463.46;
Hstar = 0.7846;
Kc = 0.1137;
tspan = 0:0.01:40;
initial_conditions = [0;0;0;0];
function dy=pair1(t,y,a,b,d,hstress,sigma,Hstar,kc,c) % local function definition
dy=zeros(4,1);
dy(1) = a*sinh((b*(sigma^d)*(1-y(2)))/((1-y(3))*(1-y(4)))) ;
dy(2) = ((hstress*(a*sinh(b*(sigma^d)*(1-y(2)))/(1-y(3))*(1-y(4))))/(sigma^d))*(1-y(2)/Hstar);
dy(3) = (kc/3)*power((1-y(3)),4);
dy(4) = c*y(1);
end
% end % <- this end ends the function myObjective so everything after this line confuses matlab
[t,y] = ode45(@(t,y)pair(t,y,a,b,d,hstress,sigma,Hstar,Kc,c),tspan,initial_conditions) %function call
plot(t,y(:,1)) % while testing OK...
diffsquare(:,i) = ((expdata(i,2) - y(1)).^2)/n;
end
chisquare = sum(diffsquare);
global c; % Avoid globals like the plague -.even if this is for an acceptable cause put that comment in to discourage use
c = (chisquare);
Then you have to make minor modifications to your script as well:
% Here you load the data-file, now your Objective-function is applyable for multiple events and cases,
% this script file is what you copy-n-modify for different data-sets, much more reproducible and leaner
% work process.
expdata = load('ru.txt');
straindata = expdata(:,2);
Timedata = expdata(:,1);
P0 = [7E-6; 5E-6; 0.5; 0.43];
lb = [1E-6; 5E-8; 0; 0.42];
ub = [3E-5; 9E-6; 1; 0.5];
P = fmincon(@(P) myObjective(P,straindata,Timedata),P0,[],[],[],[],lb,ub);
plot(t,Straindata,'+-','Linewidth',1.5,'color','red');
hold on
plot(t,y(1),'--','Linewidth',1.5,'color','blue');
xlabel('Time');
ylabel('Strain');
Haven't checked anything but these changes you should do first.
HTH
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!