Optimization of initial conditions for better Rsquared

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.

For example the 1st run before the while loop for FR is a 4095x1 double when the while loop is on the data points for FR decreases to 3640x1 double not sure why and this causes a crash bc the fit function for R squared is no longer equal lengths.

No I dont want to compare any of the runs I want to optimize the intitial conditions L,C, and A (for pendulum 1 and 2) to match the data from the Excel CFD as close as possible. For example I start with a very low initial angle (A = 1) and low damping (c= 0.001) and short length ( L=0.1) (there is an L1 and L2 for each pendulum as well as all other variables A1 A2 etc.) I want to increase or decrease these in order to create a better fit and increase my Rsquared if I guess and check them I get close curves but I would like to create a script to maniupalte these variables in order to increase Rsquared and then plot the final curve from model against excel CFD data

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])