Optimizing R squared by changing initial conditions and rerunning model
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Respuestas (0)
Ver también
Categorías
Más información sobre Computational Fluid Dynamics (CFD) en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!