clc; clear; close all;
%% 1. 风力机简化模型（2 MW）
J   = 4.5e6;  R  = 40;   ρ  = 1.225;
k1  = 0.5*ρ*pi*R^5*0.44/ω_nom^2;   % C_p(λ)≈k1·(1/β-0.022β)
ω_nom = 2*pi*1800/60;              % 额定转速 1800 rpm
P_nom = 2e6;                       % 2 MW

%% 2. RBF-PID 参数
hidden = 6;            % 隐含层节点
eta    = 0.05;         % 学习率
alpha  = 0.1;          % 动量因子
xite   = 0.05;         % PID 学习率

%% 3. 初始化
Kp=0.5; Ki=0.01; Kd=0.001;
w  =rand(hidden,1);  c =randn(hidden,3); b=ones(hidden,1);
e1=0; e2=0; e3=0; u1=0; y1=0;

%% 4. 仿真
T = 5000; ts = 0.01; ω_ref = ω_nom*ones(T,1);
ω = zeros(T,1); β = zeros(T,1);
for k=2:T
    % 风速阶跃
    v = 8 + 2*sign(k-1500) + 0.5*randn;
    
    % 风力机模型
    Cp  = 0.44*(1/β(k-1)-0.022*β(k-1));
    Pw  = 0.5*ρ*R^2*v^3*Cp;
    ω(k) = ω(k-1) + (Pw/J - k1*β(k-1)/J)*ts;
    
    % RBF-PID
    e(k) = ω_ref(k) - ω(k);
    x = [e(k); e(k)-e1; e(k)-2*e1+e2];          % 输入
    h = exp(-sum((x-c).^2,2)./(2*b.^2));        % RBF
    y_hat = w'*h;                               % 辨识输出
    dydu = sum(w.*h.*(c(:,3)-x(3))./b.^2);      % 雅可比信息
    
    % PID 增量
    du = ts*(Kp*(e(k)-e1) + Ki*e(k) + Kd*(e(k)-2*e1+e2));
    β(k) = β(k-1) + du;
    β(k) = max(0,min(β(k),30));                 % 0-30°限幅
    
    % 更新 RBF 权值
    w = w + eta*(y_hat - ω(k))*h;
    
    % PID 增益更新
    dKp = xite*e(k)*(e(k)-e1)*dydu;
    dKi = xite*e(k)*e(k)*dydu;
    dKd = xite*e(k)*(e(k)-2*e1+e2)*dydu;
    Kp = Kp + dKp; Ki = Ki + dKi; Kd = Kd + dKd;
    
    % 记录
    e2=e1; e1=e(k); e3=e(k);
end

%% 5. 结果
subplot(2,1,1); plot((1:T)*ts,ω,'b',(1:T)*ts,ω_ref,'r--'); grid;
legend('实际转速','额定转速'); title('基于RBF-PID的风力机变桨距控制');
subplot(2,1,2); plot((1:T)*ts,β,'g'); grid; ylabel('桨距角 β [°]');