Optimizing R squared by changing initial conditions and rerunning model

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

[Geoff Hayes]

Wilfredo - please clarify what you mean by loose data points the second time. Do you mean that the plot from the initial run is replaced with the data from the second run? Do you want to visually compare the first, second, third, etc. run. Or do you mean something else?

[Wilfredo Huaman]

For example the 1st run before the while loop for FR is a 4095x1 double when the while 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 to FR.

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,

Willie

[Wilfredo Huaman]

I found the issue with the decreasing loop numbers but is there a way to change those intital conditions in order to maxamize Rsquared

For example is there a way to tell matlab + or - these varaibles by 0.001 as needed to maxamize the R squared value??