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,Servo3\n');  % Header row
stateLogFile = fopen('state_log.csv', 'w');  % Open another CSV file for writing
fprintf(stateLogFile, 'Time,Phi,Theta,X_Offset,Y_Offset,Zp,Gamma\n');  % 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,%.4f\n', ...
        currentTime, dist_1, dist_2, dist_3, Phi, Theta, xof1, yof1,a11, a22, a33);
    % fprintf(stateLogFile, '%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f\n', 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;
end

Answer 1

[Aravind]

Hello @shaherbano zaidi,

To tune a PID controller using MATLAB's PID Tuner app, you need a system model. If you do not have one, you must create it using input-output data from your system. Once you have a model, you can proceed with PID tuning in the app. Here are some methods you can use:

System Identification Approach

If you lack a mathematical model for your ball-on-plate system, you can generate one using MATLAB's "System Identification" Toolbox. This tool allows you to develop a model by analyzing input-output data collected from your system over time.

Steps to Create a System Model:

• Collect Data: Run your system and record the input (control signals) and output (sensor readings) over time.
• System Identification App: Use the "System Identification" Toolbox's app to fit a model to your data. This app helps identify linear and non-linear models by processing collected data.
• Estimate Plant Parameters: Within the app, select model structures like first-order, second-order, or higher-order models to best fit your system's dynamics.
• Validate the Model: Ensure the model represents your system accurately by comparing its simulated response with actual data.

Once you have a reliable system model, you can import it into the "PID Tuner" app for effective PID parameter tuning. Here are some tutorials and documentation pages for more details on using the "System Identification" Toolbox:

• Getting Started with the System Identification Toolbox: https://in.mathworks.com/help/ident/getting-started-1.html
• Overview of System Identification: https://www.mathworks.com/help/releases/R2021b/ident/gs/about-system-identification.html
• System Identification Workflow: https://www.mathworks.com/help/releases/R2021b/ident/gs/system-identification-workflow.
• Tutorial on identifying linear models using the System Identification App: https://www.mathworks.com/help/releases/R2021b/ident/gs/identify-linear-models-using-the-gui.html
• Tutorial on identifying non-linear models using the System Identification App: https://www.mathworks.com/help/releases/R2021b/ident/gs/identify-nonlinear-black-box-models-using-the-gui.html

PID Tuner with System Identification

The "PID Tuner" app also offers a feature to estimate plant parameters directly from response data, useful when you lack a detailed system model but can obtain response data.

Using PID Tuner:

• Import Response Data: Use your system's response data to let the PID Tuner estimate the plant model by capturing how your system reacts to a step input or other test signals.
• Estimate Plant: The PID Tuner app will analyze the response data to construct a simplified plant model, which can be first-order, second-order, or customized based on the data.
• Auto-Tune PID: With the estimated model, the PID Tuner automatically calculates optimal PID parameters, offering a quick and effective tuning solution.

Here are some resources and tutorials to get started with the "PID Tuner" app:

• Overview of PID Tuner app: https://www.mathworks.com/help/releases/R2021b/control/ref/pidtuner-app.html
• System Identification for PID Control: https://www.mathworks.com/help/releases/R2021b/slcontrol/ug/system-identification-of-plant-models.html
• Tutorial to estimate plant from measured or simulated response data: https://www.mathworks.com/help/releases/R2021b/slcontrol/ug/interactively-estimate-plant-from-measured-or-simulated-response-data.html
• Tutorial to import measured response data for plant estimation: https://www.mathworks.com/help/releases/R2021b/slcontrol/ug/import-measured-response-data-for-plant-estimation.html
• Tutorial to estimate plant parameters from response data: https://www.mathworks.com/help/releases/R2021b/control/getstart/interactively-estimate-plant-parameters-from-response-data.html

Heuristic Tuning Approach

If system identification is not feasible, manual tuning using heuristic methods like the Ziegler-Nichols method can be used. This method involves trial and error but is guided by established rules.

Ziegler-Nichols Method:

• Set Initial Gains: Start with zero integral and derivative gains.
• Find Critical Gain: Gradually increase the proportional gain until consistent oscillation occurs.
• Calculate PID Parameters: Use Ziegler-Nichols tuning rules to determine PID parameters based on observed oscillation period and critical gain.

By using the System Identification Toolbox to derive a model or employing heuristic methods like Ziegler-Nichols, you can effectively tune your PID controller directly in MATLAB without needing Simulink.

I hope this answers your question.