Problems using ODE45 using fmincon command for least square optimisation problem

12 views (last 30 days)
ruban ramalingam on 13 Jul 2019
Answered: ruban ramalingam on 14 Jul 2019
Script file
function chisquare=myObjective(P)
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

Bjorn Gustavsson on 13 Jul 2019
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.
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

ruban ramalingam on 14 Jul 2019
Sir,, Whether the following codes should be in seperate file??
[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);
Actually, I was solving coupled ordinary differential equations using ODE45 and I want to do determine the optimised variables with the help of fmincon command..