基于matlab实现单摆小球运动的仿真模拟，并给出详细代码，具体要求如下——起初，单摆小球在力矩的作用下稳定在与竖直方向夹角为pi/6的位置；然后对小球施加微小扰动（如正弦扰动），请设计一个PID控制（或其他控制器）使得单摆小球稳定在与竖直方向夹角为pi/6初始位置。并绘制小球的位置（角度）随时间变化的图像和利用VideoWriter函数绘制出施加扰动前后至趋于稳定时单摆小球的模拟实物（包括固定点、单摆杆和小球三部分的模拟实物）动态图。参考思路如下： 明确单摆的动力学模型 设计一个PID控制器来控制小球的位置，将PID控制器与单摆的动力学模型相结合，得到闭环控制系统的运动方程。具体而言，我们可以将单摆的动力学模型表示为一个状态空间方程：根据PID控制器的输出，我们可以得到控制输入，将上述控制输入代入状态空间方程中，可以得到闭环控制系统的运动方程。然后可以使用MATLAB来实现上述控制系统的仿真模拟。具体而言，我们可以使用ode45函数来求解状态空间方程的数值解。 请给出详细的matlab代码！

以下是基于MATLAB实现单摆小球运动的仿真模拟的详细代码，包括扰动和PID控制：

%% 定义参数
m = 0.1;    % 小球质量
l = 0.5;    % 摆杆长度
g = 9.8;    % 重力加速度
theta0 = pi/6;    % 初始角度
omega0 = 0;    % 初始角速度
A = 0.05;    % 扰动振幅
f = 2;    % 扰动频率

%% 定义控制器参数
Kp = 50;    % 比例增益
Ki = 10;    % 积分增益
Kd = 20;    % 微分增益

%% 定义状态空间方程
A = [0 1; -g/l -0.2];    % 系统矩阵
B = [0; 1/(m*l^2)];    % 输入矩阵
C = [1 0];    % 输出矩阵
D = 0;    % 直接通道

sys = ss(A,B,C,D);    % 定义状态空间系统

%% 设计PID控制器
s = tf('s');    % 定义s变量
PID = Kp + Ki/s + Kd*s;    % 定义PID控制器
C_sys = tf(PID);    % 将PID控制器转化为系统
sys_cl = feedback(C_sys*sys,1);    % 将PID控制器与单摆系统结合成闭环控制系统

%% 求解系统的数值解
tspan = [0 20];    % 时间范围
x0 = [theta0; omega0];    % 初始状态
[t,x] = ode45(@(t,x)pendulum_sys(t,x,A,B,C,Kp,Ki,Kd,m,l,g,A,f), tspan, x0);    % 求解ODE

%% 绘制小球位置随时间变化的图像
theta = x(:,1);    % 摆角
figure;
plot(t,theta);
xlabel('时间（秒）');
ylabel('摆角（弧度）');
title('小球位置随时间变化的图像');

%% 绘制模拟实物动态图
v = VideoWriter('pendulum.avi');    % 定义视频对象
open(v);    % 打开视频对象
figure;
for i = 1:length(t)
    % 获取当前时刻的状态
    theta_i = x(i,1);
    omega_i = x(i,2);
    % 计算小球的位置
    x_i = l*sin(theta_i);
    y_i = -l*cos(theta_i);
    % 绘制模拟实物
    plot([0 x_i],[0 y_i],'k-','LineWidth',2);
    hold on;
    plot(x_i,y_i,'ro','MarkerSize',10,'MarkerFaceColor','r');
    plot([-1 1],[0 0],'k-');
    axis([-1 1 -1 1]);
    title(['t = ' num2str(t(i)) 's']);
    hold off;
    % 添加当前帧到视频对象中
    frame = getframe(gcf);
    writeVideo(v,frame);
end
close(v);

%% 定义状态空间方程的函数
function dxdt = pendulum_sys(t,x,A,B,C,Kp,Ki,Kd,m,l,g,A,f)
    % 获取当前状态
    theta = x(1);
    omega = x(2);
    % 计算控制输入
    u = -Kp*theta - Ki*trapz(t,theta) - Kd*omega;
    % 添加扰动
    theta = theta + A*sin(2*pi*f*t);
    % 计算状态变化率
    dxdt = A*[omega; -g/l*sin(theta)] + B*u;
end

注意，这里的控制器是一个简单的PID控制器，可能需要根据实际情况进行调整。同时，由于扰动的存在，控制器可能需要更长的时间才能使小球稳定在初始位置。