✅作者简介：热爱科研的Matlab仿真开发者，擅长毕业设计辅导、数学建模、数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。
🍎 往期回顾关注个人主页：Matlab科研工作室
👇 关注我领取海量matlab电子书和数学建模资料
🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信。
🔥 内容介绍
一、水箱系统控制的背景与挑战
水箱系统广泛应用于工业生产、城市供水、建筑消防等诸多领域。在工业生产中，水箱为生产流程提供稳定的水源，其水位的精确控制对生产的连续性和产品质量至关重要；城市供水系统里，水箱起到储存和调节水量的作用，以满足不同时段用户的用水需求。
然而，水箱系统的控制面临着一系列挑战：
- 系统动态特性复杂：水箱的动态特性受到多种因素影响，如管道阻力、水箱泄漏、进水和出水流量的变化等。这些因素相互作用，使得水箱水位与进出水流量之间的关系呈现出非线性和时变性，传统的基于固定参数模型的控制方法难以适应这种复杂特性。
- 外部干扰多：水箱系统容易受到外部干扰，例如用水需求的突然变化（如城市供水系统中的用水高峰）、环境温度变化对水的密度和蒸发量的影响等。这些干扰会导致水箱水位波动，增加了控制的难度，要求控制系统具备较强的抗干扰能力。
二、在线数据驱动 PID 控制原理
PID 控制基础
编辑
在线数据驱动的增强
在传统 PID 控制基础上，在线数据驱动的 PID 控制利用实时获取的系统运行数据来优化控制参数。随着水箱系统的运行，持续采集与水箱状态相关的数据，如水位、进出水流量、水温等。通过数据分析方法（如自适应算法、机器学习算法）实时调整 Kp、Ki 和 Kd。例如，当水箱受到突然的用水需求增加导致水位快速下降时，在线数据驱动机制可依据新采集的数据，迅速调整 PID 参数，使控制器能更有效地应对这种干扰，维持水箱水位稳定。这种方法能够自动适应水箱系统动态特性的变化，显著提高控制的准确性和鲁棒性。
三、离线 DD - FRIT 控制原理
DD - FRIT 控制概述
DD - FRIT 即数据驱动的前馈 - 反馈迭代学习控制（Data - Driven Feed - forward and Feedback Iterative Learning Control）。它融合了前馈控制、反馈控制和迭代学习控制的理念，借助历史数据来设计控制器。之所以称为离线，是因为该控制策略在系统实际运行前，通过对历史数据的深入分析和处理来确定控制参数。
前馈控制
前馈控制基于对系统输入（如进水流量）与输出（水箱水位）之间关系的先验知识。通过建立水箱系统的数学模型（例如基于物理原理的质量守恒方程：进水流量 - 出水流量 = 水箱水位变化率 × 水箱横截面积），根据已知的输入信号（如设定的进水流量变化规律），提前计算出控制输入，以抵消预期的干扰对输出的影响。例如，如果已知某时段用水需求会大幅增加，可通过前馈控制提前增大进水流量，以维持水位稳定。前馈控制能够在干扰影响系统输出之前就采取措施进行补偿，具有前瞻性。
反馈控制
反馈控制与 PID 控制中的反馈类似，通过测量实际输出（水箱当前水位）与期望输出（设定水位）之间的误差，对控制输入进行调整。在 DD - FRIT 中，反馈控制与前馈控制相互配合，弥补前馈控制由于模型不精确或未考虑到的干扰因素导致的控制误差。当实际水位与设定水位存在偏差时，反馈控制根据这个误差调整进水流量，使水位逐渐趋近于设定值。
迭代学习控制
迭代学习控制利用系统以往运行过程中的数据和控制经验。在每次运行结束后，根据本次运行的误差，调整下一次运行的控制输入。通过多次迭代，使系统的输出逐渐接近期望输出。例如，在水箱系统的多次运行中，每次记录水位误差，然后根据这些误差信息优化下一次运行时的控制策略，最终使水箱水位能更准确地跟踪设定值。在离线阶段，通过对大量历史数据的分析和迭代计算，确定出适用于该水箱系统的控制参数，在实际运行时直接应用这些参数进行控制。
四、两种控制策略结合的优势
在线数据驱动 PID 控制能够实时响应系统变化，对即时干扰具有良好的抑制能力；而离线 DD - FRIT 控制则利用历史数据，通过深入分析和迭代优化，提供一个相对全局最优的控制框架。将两者结合，可以发挥各自优势，使水箱系统在面对各种复杂情况时，既能快速响应实时变化，又能基于历史经验保持长期稳定运行，从而实现更高效、精确的水箱系统控制。
⛳️ 运行结果
编辑
📣 部分代码
% %僞儞僋僔僗僥儉僥僗僩梡 Simulink
% clear
% close all
% clc
%
% %慡懱僒儞僾儕儞僌帪娫
% Endtime = 600;
% Ts_sys = 0.1;
% Ts = 0.1;
%
% starttime = 3;
% endtime = 600;
%
% K = [.001 .5 .9];
%
% filename = 'Tank_test_init_sim_2017b';
% % filename = 'Tank_test_init_sim';
% open(filename);
% sim(filename);
% clear
close all
clc
Endtime = 600;
Ts_sys = 1;
starttime = 3;
endtime = Endtime/Ts_sys;
runtime = 1:endtime+1;
nu = 2;
ny = 3;
sigma = 130;
delta = 0;
Ts = 1;
rho = Ts/sigma;
mu = 0.25*(1-delta) + 0.51*delta;
t1 = -2*exp(-rho/(2*mu))*cos(sqrt(4*mu-1)/(2*mu)*rho);
t2 = exp(-rho/mu);
lr = [0.003 0.00005 0]; %[ni np nd]
knn = 6;
epochs = 3;
N_0 = endtime+1;
% K = [.001 .02 0];
% K = [0.5 .01 0];
K = [0.3 .003 0];
filename = 'Tank_test_init_sim_2017b';
open(filename);
sim(filename);
u0 = logsout.getElement('u0').Values.Data;
y0 = logsout.getElement('y0').Values.Data;
r0 = logsout.getElement('r0').Values.Data;
db = [[r0(2:endtime+1); zeros(1,1)] r0 y0 [zeros(1,1); y0(1:endtime)] [zeros(1,2).'; y0(1:endtime-1)] [zeros(1,1); u0(1:endtime)] repmat(K,endtime+1,1)];
x = zeros(endtime, 1);
r1 = zeros(endtime, 1);
K_nn = ones(endtime, knn);
wi = zeros(endtime, knn);
d = zeros(endtime, N_0);
K_old = zeros(endtime+1, 3);
K_new = zeros(endtime, 3);
grad = zeros(endtime, 3);
yr = zeros(endtime, 1);
J_ep = zeros(epochs, 1);
e = zeros(endtime, 1);
% initializing K_old and K_new
K_old(1:starttime, :) = repmat(K,starttime,1);
K_new(1:starttime, :) = repmat(K,starttime,1);
% initializing N
N(starttime) = N_0;
for ep = 1:epochs
for t = starttime:endtime
% information vector / Query
query_t = [r0(t+1); r0(t); y0(t:-1:t-ny+1); u0(t-1:-1:t-nu+1)].';
% k nearest neighbours
for j = 1:N_0
d(t,j) = sum(abs((query_t - db(j,1:nu+ny+1))./(max(db(1:N_0,1:nu+ny+1)) - min(db(1:N_0,1:nu+ny+1)))));
end
[~, sorted_index] = sort(d(t,1:N_0));
K_nn(t, :) = sorted_index(1:knn);
%% Step 3 : Computing PID Parameters
wi(t, :) = exp(-d(t,K_nn(t, :)))/sum(exp(-d(t,K_nn(t, :))));
K_old(t, :) = wi(t, :)*db(K_nn(t,:), nu+ny+2:nu+ny+4);
%% Step 4 : FRIT PID Parameters Adjustment
r1(t) = y0(t) + 1/K_old(t,2)*(u0(t) - u0(t-1) + K_old(t,1)*(y0(t) - y0(t-1)) + K_old(t,3)*(y0(t) - 2*y0(t-1) + y0(t-2)));
yr(t+1) = -t1*yr(t) - t2*yr(t-1) + (1+t1+t2)*r1(t);
e(t+1) = y0(t+1) - yr(t+1);
x0_t = u0(t) - u0(t-1) + K_old(t,1)*(y0(t) - y0(t-1)) + K_old(t,3)*(y0(t) - 2*y0(t-1) + y0(t-2));
grad(t, :) = e(t+1)*(1+t1+t2)*[-(y0(t)-y0(t-1))/K_old(t,2); x0_t/K_old(t,2)^2; -(y0(t)-2*y0(t-1)+y0(t-2))/K_old(t,2)];
K_new(t, :) = K_old(t, :) - lr.*grad(t, :);
% Updating the PID parameters
db(t, nu+ny+2:nu+ny+4) = [K_new(t, 1); K_new(t,2); 0].';
end
J_ep(ep) = 1/endtime*sum(e.^2);
end
% disp('delay start');
% pause(180);
% disp('delay end');
figure(5)
plot(1:epochs, J_ep);
title('Cost function');
xlabel('epochs');
grid on;
filename = 'Tank_test_output_sim_2017b';
open(filename);
sim(filename);
u = logsout.getElement('u').Values.Data;
y = logsout.getElement('y').Values.Data;
r = logsout.getElement('r').Values.Data;
figure(1)
subplot(2,1,1);
plot(runtime, r0, '--r', runtime, y0);
ylabel('y');
xlabel('t[steps]');
legend({'r', 'fixed pid output'}, 'Location','northwest');
axis([0 600 0 200]);
grid on;
subplot(2,1,2);
plot(runtime, u0);
ylabel('u');
xlabel('t[steps]');
grid on;
figure(2)
subplot(2,1,1);
plot(runtime, r, '--r', runtime, y);
ylabel('y');
xlabel('t[steps]');
legend({'r', 'dd frit output'}, 'Location','northwest');
axis([0 600 0 200]);
grid on;
subplot(2,1,2);
plot(runtime, u);
ylabel('u');
xlabel('t[steps]');
grid on;
figure(3)
subplot(2,1,1);
plot(runtime, r0, '--r', runtime, y0, 'm', runtime, y, 'b');
ylabel('y');
xlabel('t[steps]');
legend({'r', 'fixed pid', 'dd frit'}, 'Location','northwest');
axis([0 600 0 200]);
grid on;
subplot(2,1,2);
plot(runtime, u0, 'm', runtime, u, 'b');
ylabel('u');
xlabel('t[steps]');
grid on;
figure(4)
subplot(3,1,1);
plot(runtime, K_old(:,1));
ylabel('Kp');
xlabel('t[steps]');
grid on;
subplot(3,1,2);
plot(runtime, K_old(:,2));
ylabel('Ki');
xlabel('t[steps]');
grid on;
subplot(3,1,3);
plot(runtime, K_old(:,3));
ylabel('Kd');
xlabel('t[steps]');
grid on;
