Optimizing R squared by changing initial conditions and rerunning model
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Hello I am running a 2 pendulum simulation which involves length, initial angle, viscous damping as the initial conditions (Not a double pendulum just adding both at end) the model is run and compared to an excel data sheet and an Rsquared is calculated I basically want to change those initial conditions and rerun the model while trying to maximize the R squared to at least above 0.8.
Here is the Code I have It will run the first time but loose data points the second time not sure why? I also would like to be able to tell Matlab to add or subtract from the initital conditions as needed to maximize Rsquared.
Thanks in advanced,
Willie
%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 2D pendulum motion % %
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pendulum Initial Conditions
g = 9.81;
mass=67500; % [kg] 1st pendulum bob mass
mass2=67500; % [kg] 2nd pendulum bob mass
L=0.5; % [m] 1st pendulum rod length
L2=0.5; % [m] 2nd pendulum rod length
c = 0.001;% 0.05; % 1st Damping coefficient
c2 =0.001;% 0.1; % 2nd Damping coefficient
A1 = 1; % 1st Pendulum initial angle [deg]
A2 = 1; % 2nd Pendulum initial angle [deg]
t=90; % [s] simulation time
dt= 2/91; % Time step *****(Must Match CFD)*****
i= t/dt; % Number of iterations
%**************************************%
% % 1st Pendulum Equations of Motion % %
%**************************************%
% 1st Pendulum variable allocation
theta=zeros(round(i),1);
omega=zeros(round(i),1);
alpha=zeros(round(i),1);
Fr=zeros(round(i),1);
theta(1,:)=deg2rad(A1); % initial angular position
omega(1,:)=0; % initial angular velocity
alpha(1,:)=0; % initial angular acceleration
for n=1:i
theta(n+1,:)= theta(n,:)+omega(n,:)*dt;
omega(n+1,:)= omega(n,:)+alpha(n,:)*dt;
alpha(n+1,:)=(-g*sin(theta(n+1,:)))/L-c*omega(n+1,:);
Fr(n+1,:)= mass*g*cos(theta(n+1,:))+ mass*L*(omega(n+1,:)).^2;
end
%**************************************%
% % 2nd pendulum equations of motion % %
%**************************************%
% 2nd Pendulum variable allocation
t=0:i;
theta2=zeros(round(i),1);
omega2=zeros(round(i),1);
alpha2=zeros(round(i),1);
Fr2=zeros(round(i),1);
theta2(1,:)=deg2rad(A2); % initial angular position
omega2(1,:)=0; % initial angular velocity
alpha2(1,:)=0; % initial angular acceleration
for n=1:i
theta2(n+1,:)= theta2(n,:)+omega2(n,:)*dt;
omega2(n+1,:)= omega2(n,:)+alpha2(n,:)*dt;
alpha2(n+1,:)=(-g*sin(theta2(n+1,:)))/L2-c2*omega2(n+1,:);
Fr2(n+1,:)= mass2*g*cos(theta2(n+1,:))+ mass2*L2*(omega2(n+1,:)).^2;
end
% Combined Pendulum Force
FR = (Fr + Fr2)/2;
%R squared Value curve fitting
Fit = fitlm(TotalForce,FR);
R2 = Fit.Rsquared.Ordinary;
while R2 < 0.5
dx = 0.001;
% Pendulum Initial Conditions
g = 9.81;
mass=67500; % [kg] 1st pendulum bob mass
mass2=67500; % [kg] 2nd pendulum bob mass
L= L + dx; % [m] 1st pendulum rod length
L2=L + dx; % [m] 2nd pendulum rod length
c = c + dx;% 0.05; % 1st Damping coefficient
c2 =c2 +dx;% 0.1; % 2nd Damping coefficient
A1 = A1+dx; % 1st Pendulum initial angle [deg]
A2 = A2+dx; % 2nd Pendulum initial angle [deg]
t=80; % [s] simulation time
dt= 2/91; % Time step *****(Must Match CFD)*****
i= t/dt; % Number of iterations
%**************************************%
% % 1st Pendulum Equations of Motion % %
%**************************************%
% 1st Pendulum variable allocation
theta=zeros(round(i),1);
omega=zeros(round(i),1);
alpha=zeros(round(i),1);
Fr=zeros(round(i),1);
theta(1,:)=deg2rad(A1); % initial angular position
omega(1,:)=0; % initial angular velocity
alpha(1,:)=0; % initial angular acceleration
for n=1:i
theta(n+1,:)= theta(n,:)+omega(n,:)*dt;
omega(n+1,:)= omega(n,:)+alpha(n,:)*dt;
alpha(n+1,:)=(-g*sin(theta(n+1,:)))/L-c*omega(n+1,:);
Fr(n+1,:)= mass*g*cos(theta(n+1,:))+ mass*L*(omega(n+1,:)).^2;
end
%**************************************%
% % 2nd pendulum equations of motion % %
%**************************************%
% 2nd Pendulum variable allocation
t=0:i;
theta2=zeros(round(i),1);
omega2=zeros(round(i),1);
alpha2=zeros(round(i),1);
Fr2=zeros(round(i),1);
theta2(1,:)=deg2rad(A2); % initial angular position
omega2(1,:)=0; % initial angular velocity
alpha2(1,:)=0; % initial angular acceleration
for n=1:i
theta2(n+1,:)= theta2(n,:)+omega2(n,:)*dt;
omega2(n+1,:)= omega2(n,:)+alpha2(n,:)*dt;
alpha2(n+1,:)=(-g*sin(theta2(n+1,:)))/L2-c2*omega2(n+1,:);
Fr2(n+1,:)= mass2*g*cos(theta2(n+1,:))+ mass2*L2*(omega2(n+1,:)).^2;
end
% Combined Pendulum Force
FR = (Fr + Fr2)/2;
%R squared Value curve fitting
Fit = fitlm(TotalForce,FR);
R2 = Fit.Rsquared.Ordinary;
if c && c2 > 0.9
c = 0.9;
c2 = 0.9;
elseif A1 && A2 > 50
break
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Plotting Pendulum & CFD together % %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
plot(t*dt,Fr,t*dt,Fr2)
legend('Pendulum 1','Pendulum 2')
xlim([5 90])
figure(3)
plot(t*dt,FR,time,TotalForce,'LineWidth',2)%,t*dt,Fr,t*dt,Fr2,)
legend('Resultant Pendulum','CFD')
%legend('Pendulum 1','Pendulum 2','Resultant Pendulum','CFD','location','bestoutside')
xlim([5 90])
3 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Computational Fluid Dynamics (CFD) finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!