How to make an adaptive PID controller?

The PID program I write currently works fine with the DC motor of our robot. Considering that the friction efficient and mass may vary(we need to load/unload weights to the robot), I am thinking about improving the system by using an adaptive PID controller. Does anyone have some sample code/algorithm I can refer to? I am new to this type of PID so simple tutorials are welcomed too.

Thanks in advance!

Answer 1

[Birdman]

Check the following links:

https://www.mathworks.com/matlabcentral/fileexchange/44416-simple-adaptive-control-example

https://www.mathworks.com/matlabcentral/fileexchange/18061-result-adaptive-pid-control-improves-laserscanner-positioning

Answer 2

[Raydeen]

%% =====================================================================

%% MPHS Adaptive PID + NARX Feed-forward – Three Pre-heaters in Parallel

%% Tenova South Africa (Pty) Ltd – Aksu Kazchrome Furnace 64 MPHS

%% Author: M365 Copilot for Ramjee Raydeen (TENOVA)

%% Date: 2025-12-09 (Corrected 2024-05-16)

%%

%% Description

%% Extends Modes 2/3 with:

%% - Adaptive PID (IMC-style retuning) using online RLS identification

%% - NARX feed-forward prediction to pre-empt disturbances

%% Compares three variants:

%% (1) Baseline PIDF

%% (2) Adaptive PIDF (RLS-IMC)

%% (3) Adaptive PIDF + NARX feed-forward

%%

%% Requirements

%% Control System Toolbox; Deep Learning Toolbox (for narxnet). If narxnet

%% is unavailable, a simple linear predictor is used as fallback.

%%

%% Output

%% - KPI bar charts per vessel comparing OS, Ts, IAE across variants

%% - Figures saved into ./outputs/adaptive_overview.png and adaptive_kpis.png

%% =====================================================================

clear; close all; clc; rng(44);

%% ======= Settings =======

Tsim = 1800; dt = 0.5; t = 0:dt:Tsim;

num_steps = numel(t);

mode = 2; % set 2 for Level-4 temp or 3 for Off-gas temp

% FIX: Restore the intended setpoints arrays

Tk_set = [250, 250, 250]; % Level-4 temp setpoints [°C]

Tog_set= [ 90, 90, 90]; % Off-gas temp setpoints [°C]

P_band = struct('low',0.95,'high',1.05);

PH(1) = struct('K',2.4,'tau',90,'theta',10,'name','PH1');

PH(2) = struct('K',2.1,'tau',75,'theta',12,'name','PH2');

PH(3) = struct('K',2.6,'tau',80,'theta', 9,'name','PH3');

OG(1) = struct('a',0.18,'b',0.09,'c',35);

OG(2) = struct('a',0.16,'b',0.11,'c',35);

OG(3) = struct('a',0.17,'b',0.10,'c',35);

u_min = 0; u_max = 100; du_max = 8;

manifold = struct('tauP',25,'Kp',0.015,'c',0.006,'P0',1.00);

Power_min = 4; Power_max = 18; % MW

GcP = pidtune(tf(1, [manifold.tauP 1]), 'PI');

measNoiseStd_T = 0.5; measNoiseStd_OG = 0.4;

dist = zeros(3, numel(t));

dist(:, t>=600 & t<900) = dist(:, t>=600 & t<900) + repmat([10; -7; 5],1,sum(t>=600 & t<900));

dist(:, t>=1200) = dist(:, t>=1200) + repmat([-8; 9; -6],1,sum(t>=1200));

plants = cell(1,3); C0 = cell(1,3);

for i=1:3

plants{i} = tf(PH(i).K, [PH(i).tau 1], 'InputDelay', PH(i).theta);

opts = pidtuneOptions('DesignFocus','balanced');

C0{i} = pidtune(plants{i}, 'PIDF', 1/PH(i).tau, opts);

C0{i}.Kd = 0.5*C0{i}.Kd; C0{i}.Tf = max(4*C0{i}.Tf, C0{i}.Tf);

end

%% ======= NARX preparation (offline synthetic training; replace with plant data) =======

useNarx = (exist('narxnet','file') == 2);

Kff = [0.4, 0.4, 0.4]; % feed-forward blend

narxDelays = 1:4; narxHidden = 8;

[Utrain, Ytrain] = synthNarxTraining(t);

if useNarx

nets = cell(1,3);

for i=1:3

% FIX: Specify number of inputs (4 channels) explicitly in narxnet definition

nets{i} = narxnet(narxDelays, narxDelays, narxHidden, 'open', 4);

nets{i}.trainParam.showWindow = false;

% Prepare training data: Utrain is (4 x T) data, Ytrain is (1 x T) data

[Xs,Xi,Ai,Ts] = preparets(nets{i}, num2cell(Utrain, 1), {}, num2cell(Ytrain, 1));

nets{i} = train(nets{i}, Xs, Ts, Xi, Ai);

% Convert to closed loop for prediction in the simulation loop

nets{i} = closeloop(nets{i});

end

% Pre-allocate state variables for sequential simulation

Xi_state = cell(1,3);

Ai_state = cell(1,3);

else

nets = []; % fallback will be used in simulation

end

%% ======= Run variants =======

variants = {'baseline','rls','rls_narx'}; Nv = numel(variants);

metricsOS = zeros(3,Nv); metricsTs = zeros(3,Nv); metricsIAE = zeros(3,Nv);

% Ensure 'outputs' directory exists once at the start

outputsDir = fullfile(pwd, 'outputs');

if ~exist(outputsDir, 'dir')

mkdir(outputsDir);

fprintf('Created directory: %s\n', outputsDir);

end

for v=1:Nv

fprintf('Starting simulation for variant: %s\n', variants{v});

% Reset all persistent states in utility functions before each variant run

clear discretePI pidf_discrete rlsUpdate

% Pre-allocate logs for this run

u = zeros(3,length(t)); T4 = zeros(3,length(t)); OGT = zeros(3,length(t));

Pman = zeros(1,length(t)); Power= zeros(1,length(t)); flow = zeros(3,length(t));

Pman(1) = manifold.P0; Power(1) = 8; T4(:,1) = [120; 115; 130]; u(:,1) = [20; 20; 20];

% Copy initial controllers (C0) to working controllers (C)

C = cellfun(@(x) x, C0, 'UniformOutput', false);

% Setup input history buffer for dead time implementation (FIXED)

max_delay_steps = ceil(max([PH(:).theta]) / dt);

u_hist_buffer = zeros(3, max_delay_steps + 1);

for i = 1:3

delay_steps = ceil(PH(i).theta / dt);

u_hist_buffer(i, (end-delay_steps):end) = u(i,1);

end

% Initialize NARX states for this variant simulation

if useNarx && strcmp(variants{v},'rls_narx')

for i = 1:3

% Use a segment of the training data to get valid initial states

U_init_segment = Utrain(:, 1:max(narxDelays)+1);

Y_init_segment = Ytrain(:, 1:max(narxDelays)+1);

[~, Xi_state{i}, Ai_state{i}, ~] = preparets(nets{i}, num2cell(U_init_segment, 1), {}, num2cell(Y_init_segment, 1));

end

end

for k=2:length(t)

% --- Supervisory cascade (Manifold Pressure) ---

P_err = 0;

if Pman(k-1) < P_band.low, P_err = (P_band.low - Pman(k-1)); end

if Pman(k-1) > P_band.high, P_err = (P_band.high - Pman(k-1)); end

dPower_cmd = discretePI(P_err, dt, GcP.Kp, GcP.Ki, 'Pman_Ctrl');

Power(k) = saturate(Power(k-1) + dPower_cmd, Power_min, Power_max);

% Manifold dynamics calculation

sumU_prev = sum(u(:,k-1));

dP_dt = (1/manifold.tauP) * (manifold.Kp * Power(k-1) - (Pman(k-1) - manifold.P0) - manifold.c * sumU_prev);

Pman(k) = Pman(k-1) + dP_dt * dt;

for i=1:3

flow_i = max(0, min(1, 0.6*(Power(k)/Power_max) * (1.05/Pman(k)) * (u(i,k-1)/100)));

flow(i,k) = flow_i;

Tog_i = OG(i).a*T4(i,k-1) + OG(i).b*(100*flow_i) + OG(i).c + measNoiseStd_OG*randn();

OGT(i,k) = Tog_i;

if mode==2, y_meas = T4(i,k-1) + measNoiseStd_T*randn(); r = Tk_set(i);

else, y_meas = Tog_i; r = Tog_set(i); end

e = r - y_meas;

% === Adaptation ===

if ~strcmp(variants{v},'baseline')

phi = [u(i,k-1); y_meas];

theta = rlsUpdate(phi, y_meas, dt, sprintf('RLS%d',i));

Khat = max(0.5, theta(1)); tauhat = max(30, 1/abs(theta(2)));

[Kp_i,Ki_i,Kd_i,Tf_i] = imcPID(Khat, tauhat, 0.15*tauhat);

C{i}.Kp = Kp_i; C{i}.Ki = Ki_i; C{i}.Kd = 0.5*Kd_i; C{i}.Tf = max(Tf_i, C{i}.Tf);

end

% === Feed-forward (NARX) ===

u_ff = 0;

if strcmp(variants{v},'rls_narx')

if useNarx

% Input vector Uff needs to be a cell array of 1 time step

% Uff_vec is (4x1)

Uff_vec = [Power(k)/Power_max; Pman(k); u(i,k-1)/100; y_meas];

Uff_cell = { Uff_vec };

% Pass current inputs and previous states, retrieve new states

[y_pred_cell, Xi_state{i}, Ai_state{i}] = nets{i}(Uff_cell, Xi_state{i}, Ai_state{i});

y_pred = cell2mat(y_pred_cell); % Extract numeric value

else

y_pred = 0.2*(Power(k)/Power_max) + 0.15*Pman(k) + 0.4*y_meas + 50; % simple fallback

end

u_ff = Kff(i) * (r - y_pred);

if Pman(k) > P_band.high || Pman(k) < P_band.low, u_ff = 0.5*u_ff; end

end

du_cmd = pidf_discrete(C{i}, e, dt, sprintf('C%d', i)) + u_ff;

u_cmd_raw = u(i,k-1) + limitRate(du_cmd, du_max, dt);

u(i,k) = saturate(u_cmd_raw, u_min, u_max);

delay_steps = ceil(PH(i).theta / dt);

u_delayed = u_hist_buffer(i, end - delay_steps);

T4(i,k) = fopdt_discrete(T4(i,k-1), u_delayed, PH(i).K, PH(i).tau, dist(i,k), dt);

end

u_hist_buffer(:, 1:end-1) = u_hist_buffer(:, 2:end);

u_hist_buffer(:, end) = u(:,k);

end

% KPIs per vessel

for i=1:3

if mode==2, y = T4(i,:); r = Tk_set(i); else, y = OGT(i,:); r = Tog_set(i); end

info = stepinfo(y, t, r);

metricsOS(i,v) = info.Overshoot; metricsTs(i,v) = info.SettlingTime;

metricsIAE(i,v)= trapz(t, abs(r - y));

end

% Plot overview for this variant

tile = tiledlayout(2,2,'TileSpacing','compact');

sgtitle(tile, sprintf('Adaptive Variant: %s | Mode %d', variants{v}, mode));

for i=1:3

nexttile;

if mode==2

y = T4(i,:); r = Tk_set(i); yl='Level-4 Temp [\circC]';

else

y = OGT(i,:); r = Tog_set(i); yl='Off-gas Temp [\circC]';

end

plot(t,y,'b','LineWidth',1.2); hold on; yline(r,'--r','Setpoint'); grid on; xlabel('Time [s]'); ylabel(yl);

title(sprintf('%s', PH(i).name)); info = stepinfo(y,t,r);

yline(r*(1 + info.Overshoot/100), ':k', sprintf('OS=%.1f%%', info.Overshoot));

xline(info.SettlingTime, ':m', sprintf('Ts=%.0fs', info.SettlingTime));

end

nexttile; plot(t,Pman,'k',t,Power,'--g','LineWidth',1.2); grid on; legend('P_{man}','Power');

yline(P_band.low, ':r','P low'); yline(P_band.high, ':r','P high'); xlabel('Time [s]'); title('Manifold & Power');

exportgraphics(tile, fullfile(outputsDir, sprintf('adaptive_overview_%s.png', variants{v})), 'Resolution', 180);

end

%% ======= KPI bar charts =======

labels = {'Baseline','RLS-IMC','RLS-IMC+NARX'};

figure('Position',[100 100 1000 400]);

tiledlayout(1,3);

for i=1:3

nexttile; bar([metricsOS(i,:); metricsTs(i,:); metricsIAE(i,:)]');

title(sprintf('KPIs – %s', PH(i).name)); grid on; xticklabels(labels);

legend({'Overshoot [%]','Ts [s]','IAE'}, 'Location','bestoutside');

end

exportgraphics(gcf, fullfile(outputsDir,'adaptive_kpis.png'), 'Resolution', 180);

%% ======= Local utilities (Corrected functions) =======

function y = saturate(x, xmin, xmax), y = min(max(x,xmin),xmax); end

function du = limitRate(du_cmd, du_rate, dt), du = max(min(du_cmd, du_rate*dt), -du_rate*dt); end

function [uTrainData, yTrainData] = synthNarxTraining(t)

% Generates synthetic training data as a 4xT matrix

T = length(t);

% Input 1: Power scaled, Input 2: Pman scaled, Input 3: u scaled, Input 4: y_meas scaled (dummy data)

uTrainData = rand(4, T);

% Target Output 1: Simple linear combination with some dynamics

yTrainData = filter(1, [1 -0.5 0.2], sum(uTrainData, 1)) + 0.1*randn(1, T);

end

function [Kp, Ki, Kd, Tf] = imcPID(K, tau, lambda)

Kp = (tau) / (K * lambda);

Ki = Kp / tau;

Kd = 0;

Tf = 0;

end

function xk = fopdt_discrete(xkm1, u_delayed, K, tau, d, dt)

dx_dt = (-xkm1 + K * (u_delayed/100)) / tau;

xk = xkm1 + dx_dt * dt + 0.1 * d * dt;

end

function du = discretePI(e, dt, Kp, Ki, ID)

persistent I_state_map

if isempty(I_state_map), I_state_map = containers.Map('KeyType','char', 'ValueType','double'); end

if ~isKey(I_state_map, ID), I_state_map(ID) = 0; end

I_state_map(ID) = I_state_map(ID) + Ki * e * dt;

du = Kp * e + I_state_map(ID);

end

function du_cmd = pidf_discrete(C, e, dt, ID_str)

persistent I_state_map D_state_map

if isempty(I_state_map), I_state_map = containers.Map('KeyType','char', 'ValueType','double'); end

if isempty(D_state_map), D_state_map = containers.Map('KeyType','char', 'ValueType','double'); end

if ~isKey(I_state_map, ID_str), I_state_map(ID_str) = 0; end

if ~isKey(D_state_map, ID_str), D_state_map(ID_str) = 0; end

Kp = C.Kp; Ki = C.Ki; Kd = C.Kd; Tf = C.Tf;

I_state_map(ID_str) = I_state_map(ID_str) + Ki * e * dt;

D_new = (Tf/(Tf+dt))*D_state_map(ID_str) + (Kd*dt/(Tf+dt))*e/dt;

D_state_map(ID_str) = D_new;

du_cmd = Kp * e + I_state_map(ID_str) + D_state_map(ID_str);

end

function theta_hat = rlsUpdate(phi, y_meas, dt, ID_str)

persistent P_map Theta_map Lambda

if isempty(P_map), P_map = containers.Map('KeyType','char', 'ValueType','any'); end

if isempty(Theta_map), Theta_map = containers.Map('KeyType','char', 'ValueType','any'); end

if isempty(Lambda), Lambda = 0.99; end

if ~isKey(P_map, ID_str)

P_map(ID_str) = eye(length(phi)) * 100;

Theta_map(ID_str) = [0.1; 0.1];

end

P = P_map(ID_str);

Theta = Theta_map(ID_str);

K_gain = P * phi / (Lambda + phi' * P * phi);

e_pred = y_meas - phi' * Theta;

Theta = Theta + K_gain * e_pred;

P = (P - K_gain * phi' * P) / Lambda;

P_map(ID_str) = P;

Theta_map(ID_str) = Theta;

theta_hat = Theta;

end