How to make matlab and simulink work simultaneously.

Hi, im working on a ball on a plate system, i dont have a system model and i want to auto tune my PID controller on the basis of the data received from sensors. I want to know how to auto tune PID in matlab or if theres a way to connect sensor data with simulink in real time for auto tuning purpose. Moreover, i find it difficult to perfrom image processing to gather data in simulink. My code is shared below for further assistance.
%% Initialize Environment
clear all;
clc;
close all;
%% Initialize Data Logging
logFile = fopen(‘data_log.csv’, ‘w’); % Open a CSV file for writing
fprintf(logFile, ‘Time,Disturbance1,Disturbance2,Disturbance3,RollAngle,PitchAngle,CentroidX,CentroidY,Servo1,Servo2,Servo3n’); % Header row
stateLogFile = fopen(‘state_log.csv’, ‘w’); % Open another CSV file for writing
fprintf(stateLogFile, ‘Time,Phi,Theta,X_Offset,Y_Offset,Zp,Gamman’); % Header row

%% Initialize Image Processing
% Load the camera parameters
load(‘test.mat’); % Ensure the file ‘cameraParams.mat’ is in your working directory
% Create a webcam object
cap = webcam(‘Adesso CyberTrack H3’);
dist_1 = 0;
dist_2=dist_1
dist_3=dist_2
% Adjust camera settings if supported
try
% Set exposure time within valid range
exposureTime = -6; % Adjust as needed within the valid range
cap.ExposureMode = ‘manual’;
cap.Exposure = exposureTime;
catch ME
disp(‘Warning: Unable to set exposure time for the camera.’);
disp(ME.message);
end

%% Anomaly Variables
% Parameters for angle adjustments
p_a = 0;
p_r = 0;

%% Blob Object
% Blob analysis object for ball detection
% Blob analysis object for ball detection using vision.BlobAnalysis
blobAnalyzer = vision.BlobAnalysis(…
‘MinimumBlobArea’, 200, … % Adjust to filter small noise
‘MaximumBlobArea’, 50000, … % Adjust to ensure the ball is detected
‘BoundingBoxOutputPort’, true, …
‘CentroidOutputPort’, true, …
‘AreaOutputPort’, true);

%% Create Figure and Axes
hFig = figure(‘DoubleBuffer’, ‘on’);
hAxes = axes(‘Parent’, hFig);
s = 1;
%% Main Loop

%% Measure Sampling Time
%% Servo Motor Control Initial Position
% Port at which your Arduino is connected
port = ‘COM5’;
% Model of your Arduino board
board = ‘Uno’;
% Creating Arduino object with servo library
arduino_board = arduino(port, board, ‘Libraries’, ‘Servo’);
% Creating servo motor object
servo_motor = servo(arduino_board, ‘D9′,’MinPulseDuration’, 500e-6, ‘MaxPulseDuration’, 2500e-6);
writePosition(servo_motor, 0.333); % 45 degrees in as per calculations and in radians = 0.333
servo_motor1 = servo(arduino_board, ‘D10′,’MinPulseDuration’, 500e-6, ‘MaxPulseDuration’, 2500e-6);
writePosition(servo_motor1, 0.333);
servo_motor2 = servo(arduino_board, ‘D11′,’MinPulseDuration’, 500e-6, ‘MaxPulseDuration’, 2500e-6);
writePosition(servo_motor2, 0.333);

A= 45;
B= 45;
C = 45;
disp(‘done here1’);
h = 1;
for i = 1:2
%% Image Processing Started
start = tic;
[Zp,alpha] = forwardkinfunc(A,B,C);
% Read frame from webcam
frame = snapshot(cap);

% Undistort the frame using the camera parameters
undistortedFrame = undistortImage(frame, cameraParams2);

% Resize the undistorted frame to a lower resolution (e.g., half the original size)
undistortedFrame = imresize(undistortedFrame, 0.2); % Reduce resolution to 50%
%% Calculate the center of the image
[frameHeight, frameWidth, ~] = size(undistortedFrame);
centerX = frameWidth / 2;
centerY = frameHeight / 2;

%% Circle (Dot) Detection
% Convert to grayscale and adjust contrast
I_gray = rgb2gray(undistortedFrame);
I_contrast = imadjust(I_gray);

% Detect circles in the image (dots)
[Centers, Radii] = imfindcircles(I_contrast, [1 5], ‘ObjectPolarity’, ‘dark’, ‘Sensitivity’, 0.9);

% Check if circles are detected
if isempty(Centers)
disp(‘No circles detected’);
else
disp([‘Detected ‘, num2str(size(Centers, 1)), ‘ circles’]);
end

% Display the undistorted frame and detected circles
imshow(undistortedFrame, ‘Parent’, hAxes);
hold(hAxes, ‘on’);
viscircles(hAxes, Centers, Radii, ‘EdgeColor’, ‘r’);
% Plot a magenta dot at the center of the image
plot(hAxes, centerX, centerY, ‘m.’, ‘MarkerSize’, 20); % Magenta dot

% Check if at least 2 circles are detected for line fitting
if size(Centers, 1) >= 2
% Fit lines to X and Y coordinates separately
x = Centers(:, 1);
yof = Centers(:, 2);
pX = polyfit(yof, x, 1); % Line fit for X as a function of Y
pY = polyfit(x, yof, 1); % Line fit for Y as a function of X

% Display fitted line coefficients
disp([‘Line fit for X (Y): slope = ‘, num2str(pX(1))]);
disp([‘Line fit for Y (X): slope = ‘, num2str(pY(1))]);

% Calculate the slopes
slopeX = pX(1);
slopeY = pY(1);

% Calculate roll and pitch angles from the slopes
Theta = -atan2(slopeX, 1); % Tilt along X-axis (roll)
Phi = atan2(slopeY, 1); % Tilt along Y-axis (pitch)

% Display calculated angles
disp([‘Roll Angle: ‘, num2str(Theta)]);
disp([‘Pitch Angle: ‘, num2str(Phi)]);

% Draw the X and Y axis lines
plot(hAxes, [0, size(I_contrast, 2)], [mean(yof), mean(yof)], ‘g–‘, ‘LineWidth’, 2); % X-axis line
plot(hAxes, [mean(x), mean(x)], [0, size(I_contrast, 1)], ‘y–‘, ‘LineWidth’, 2); % Y-axis line

% Draw the slope lines
yFitLineX = polyval(pX, [0, size(I_contrast, 1)]);
plot(hAxes, yFitLineX, [0, size(I_contrast, 1)], ‘b-‘, ‘LineWidth’, 2); % Slope line for X

xFitLineY = polyval(pY, [0, size(I_contrast, 2)]);
plot(hAxes, [0, size(I_contrast, 2)], xFitLineY, ‘r-‘, ‘LineWidth’, 2); % Slope line for Y

% Display the calculated angles
title(hAxes, sprintf(‘Roll: %.2f°, Pitch: %.2f°’, Theta, Phi));
else
title(hAxes, ‘Not enough circles detected’);
end

%% Ball Detection Using vision.BlobAnalysis
% Check if any blobs are detected
%% Ball Detection Using vision.BlobAnalysis
hsvFrame = rgb2hsv(undistortedFrame);

% Define thresholds for channel 1 based on histogram settings
channel1Min = 0.845;
channel1Max = 0.081;

% Define thresholds for channel 2 based on histogram settings
channel2Min = 0.133;
channel2Max = 1.000;

% Define thresholds for channel 3 based on histogram settings
channel3Min = 0.519;
channel3Max = 1.000;

% Create a binary mask based on chosen histogram thresholds
sliderBW = ((hsvFrame(:,:,1) >= channel1Min) | (hsvFrame(:,:,1) <= channel1Max)) & …
(hsvFrame(:,:,2) >= channel2Min) & (hsvFrame(:,:,2) <= channel2Max) & …
(hsvFrame(:,:,3) >= channel3Min) & (hsvFrame(:,:,3) <= channel3Max);
BW = sliderBW;
BW = imfill(BW, ‘holes’); % Fill any small holes
BW = bwareaopen(BW, 200); % Remove small noise (adjust 200 as needed)
% Apply morphological operations to clean up the mask
BW = imopen(BW, strel(‘disk’, 5));

BW = imclose(BW, strel(‘disk’, 20));

% Perform blob analysis
[area,centroid,bbox] = step(blobAnalyzer,BW);
if ~isempty(centroid)
x_p = centroid(1);
y_p = centroid(2);

% Draw the bounding box on the undistorted frame
rectangle(‘Position’, bbox, ‘EdgeColor’, ‘green’, ‘LineWidth’, 2);
plot(hAxes, x_p, y_p, ‘ro’, ‘MarkerSize’, 10);
pix_cm = 32/frameWidth;
xc = x_p*pix_cm;
yc = y_p*pix_cm
[xof1, yof1] = camera_base(xc,yc)

end
elapsed_time = toc(start); % Compute elapsed time dynamically
Ts = getElapsedTime();% Use elapsed time for PID calculations
tic;
Kp_x = 10; Ki_x = 0; Kd_x = 0
Kp_y = 10; Ki_y = 0; Kd_y = 0;
Kp_phi = 80; Ki_phi = 0; Kd_phi = 0.5;
Kp_theta =80; Ki_theta =0; Kd_theta = 0.5;

% Setpoints (target values)
setpoint_x = -0.1;
setpoint_y = -5.5;
setpoint_phi = 0.00;
setpoint_theta = 0.00;

previous_errorx = 0;
integral_errorx=0
derivative_errorx = 0

previous_errory = 0;
integral_errory=0
derivative_errory = 0

previous_errorphi = 0;
integral_errorphi=0
derivative_errorphi = 0

previous_errortheta = 0;
integral_errortheta=0
derivative_errortheta = 0
prev_filtest = 0;
filte_est =0;
prev_filtesty = 0;
filte_esty =0;
prev_filtestp = 0;
filte_estp =0;
prev_filtestth = 0;
filte_estth =0;
o = 0.05;
K_xp = 4/100;
K_yp= 4.8/100;
K_xr = 0.63/100;
K_yr = 3.05/100;

error_x= setpoint_x – xof1;
integral_errorx = integral_errorx + error_x * Ts;
derivative_errorx = (error_x – previous_errorx);
filte_est = (o*prev_filtest)+(1-o)*derivative_errorx
derivative = filte_est/Ts;
prev_filtest = filte_est;
control_signalx = Kp_x * error_x + Ki_x * integral_errorx + Kd_x * derivative;
previous_errorx = error_x;

% Cross-coupling effect on phi due to x_error
theta_coupling_effect = K_xr * error_x;

error_y= setpoint_y – yof1;
integral_errory = integral_errory + error_y * Ts;
derivative_errory = (error_y – previous_errory)
filte_esty = (o*prev_filtesty)+(1-o)*derivative_errory
derivativey = filte_esty/Ts;
prev_filtesty = filte_esty;
control_signaly = Kp_y * error_y + Ki_y * integral_errory + Kd_y * derivativey;
previous_errory = error_y;

% Cross-coupling effect on theta due to y_error
phi_coupling_effect = K_yp * error_y;

% Step 5: Compute final tilt angles
theta_desired =(K_yr * error_y) + phi_coupling_effect; % Roll (tilt about X-axis)
phi_desired = (K_xp * error_x) + theta_coupling_effect; % Pitch (tilt about Y-axis)

error_phi= phi_desired – Phi;
integral_errorphi = integral_errorphi + error_phi * Ts;
derivative_errorphi = (error_phi – previous_errorphi)
filte_estp = (o*prev_filtestp)+(1-o)*derivative_errorphi
derivativep = filte_estp/Ts;
prev_filtestp = filte_estp;
control_signalphi = Kp_phi * error_phi + Ki_phi * integral_errorphi + Kd_phi * derivativep;
previous_errorphi = error_phi;
Phi1 = (Phi + control_signalphi * Ts); % Update tilt
Phi1 = max(-0.6, min(Phi1, 0.6));
pidTuner(Phi1)

error_theta= theta_desired – Theta;
integral_errortheta = integral_errortheta + error_theta * Ts;
derivative_errortheta = (error_theta – previous_errortheta)
filte_estth = (o*prev_filtestth)+(1-o)*derivative_errortheta;
derivativeth = filte_estth/Ts;
prev_filtestth = filte_estth;
control_signaltheta = Kp_theta * error_theta + Ki_theta * integral_errortheta + Kd_theta * derivativeth;
previous_errortheta = error_theta;
Theta1 = (Theta + control_signaltheta * Ts); % Update tilt
Theta1 = max(-0.6, min(Theta1, 0.6));
xof1 = xof1+(control_signalx*Ts);
yof1 = yof1+(control_signaly*Ts)

disp(‘done here 3’)
xo= 82*sin(Phi);
yo =-82*sin(Theta);

[a11f,a22f,a33f] = inversekin(xof1,yof1,Phi1,Theta1,Zp,alpha)
% [a11f,a22f,a33f] = inversekin(xof1,yof1,phi_control,theta_control,Zp,gamma)
a11 = a11f;
a22 = a22f;
a33 = a33f;
disp(‘done here 4’)
% Clamp the angles to the range 40° to 50°
a11 = max(42, min(a11, 48));
a22 = max(42, min(a22, 48));
a33 = max(42, min(a33, 48));
currentTime = datestr(now, ‘HH:MM:SS.FFF’); % Get current time
fprintf(logFile, ‘%s,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.4f,%.4f,%.4fn’, …
currentTime, dist_1, dist_2, dist_3, Phi, Theta, xof1, yof1,a11, a22, a33);
% fprintf(stateLogFile, ‘%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3fn’, Ts, Phi, Theta, xof1, yof1, Zp, gamma);

%% Servo Motor Code
% Clamp the angles to the range 40° to 50°
a = normalize(a11);
b = normalize(a22);
c = normalize(a33);

writePosition(servo_motor,a)
writePosition(servo_motor1,b)
writePosition(servo_motor2,c)
stop_it = toc(start);
A = a11;
B = a22;
C= a33;
stop = toc(start)
disp(‘done here 5’);

if (xof1 <=0.5 && xof1>=-0.5)&&(yof1 <= 1.5 &&yof1>=-1.5)
% if (xof1==0 && yof1==0)
disp(‘im done’)
writePosition(servo_motor,0.3300)
writePosition(servo_motor1,0.3300)
writePosition(servo_motor2,0.3300)
break;
else
disp(‘continued’)
h = h+1
end

end

% Clean up
fclose(logFile); % Close the log file
clear arduino_board servo_motor servo_motor1 servo_motor2 s
clear cap;
release(blobAnalyzer);
close(hFig); % Close the figure when done
disp(‘System is Stable’);

function norm_position = normalize(angle_deg)
angle_deg = angle_deg/1.5;
angle_deg = 90 – angle_deg;
angle_deg = angle_deg – (angle_deg/100);
norm_position = angle_deg/180;
endHi, im working on a ball on a plate system, i dont have a system model and i want to auto tune my PID controller on the basis of the data received from sensors. I want to know how to auto tune PID in matlab or if theres a way to connect sensor data with simulink in real time for auto tuning purpose. Moreover, i find it difficult to perfrom image processing to gather data in simulink. My code is shared below for further assistance.
%% Initialize Environment
clear all;
clc;
close all;
%% Initialize Data Logging
logFile = fopen(‘data_log.csv’, ‘w’); % Open a CSV file for writing
fprintf(logFile, ‘Time,Disturbance1,Disturbance2,Disturbance3,RollAngle,PitchAngle,CentroidX,CentroidY,Servo1,Servo2,Servo3n’); % Header row
stateLogFile = fopen(‘state_log.csv’, ‘w’); % Open another CSV file for writing
fprintf(stateLogFile, ‘Time,Phi,Theta,X_Offset,Y_Offset,Zp,Gamman’); % Header row