参照以上代码，基于matlab实现单摆小球运动的仿真模拟，并给出详细代码，具体要求如下——起初，单摆小球在力矩的作用下稳定在与竖直方向夹角为pi/6的位置；然后对小球施加微小扰动（如正弦扰动），请设计一个PID控制（或其他控制器）使得单摆小球稳定在与竖直方向夹角为pi/6初始位置。并绘制小球的位置（角度）随时间变化的图像和利用VideoWriter函数绘制出施加扰动前后至趋于稳定时单摆小球的模拟实物（包括固定点、单摆杆和小球三部分的模拟实物）动态图。

下面是基于MATLAB的单摆小球运动的仿真模拟的代码。其中，我们使用PID控制器来控制单摆小球的运动状态，使其从扰动后的状态回到初始状态。

clear all
clc

global m l g
m = 1; l = 1; g = 9.8;

% 初始状态
theta0 = pi/6; % 初始角度
omega0 = 0; % 初始角速度

% 扰动信号（正弦信号）
t = linspace(0, 10, 1000);
disturbance = 0.1*sin(5*t);

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

% 模拟时间和时间步长
tspan = [0 10];
dt = 0.01;

% 初始状态向量
y0 = [theta0; omega0];

% 求解微分方程组
[t, y] = ode45(@ode_pendulum, tspan, y0);

% PID控制器
theta_desired = theta0; % 目标角度
theta_error = zeros(size(t)); % 角度误差
theta_integral = 0; % 角度误差积分
theta_derivative = 0; % 角度误差微分
tau = zeros(size(t)); % 控制力矩

% 循环计算PID控制器的输出
for i = 1:length(t)
    % 计算角度误差
    theta_error(i) = theta_desired - y(i, 1);
    
    % 计算角度误差积分
    theta_integral = theta_integral + theta_error(i) * dt;
    
    % 计算角度误差微分
    if i > 1
        theta_derivative = (theta_error(i) - theta_error(i-1)) / dt;
    end
    
    % 计算控制力矩
    tau(i) = Kp * theta_error(i) + Ki * theta_integral + Kd * theta_derivative;
end

% 绘制角度随时间变化的曲线
figure(1)
plot(t, y(:, 1), 'r', t, theta_desired*ones(size(t)), 'b')
xlabel('Time')
ylabel('Angle')
title('Angle vs. Time')
legend('Angle', 'Desired Angle')

% 绘制控制力矩随时间变化的曲线
figure(2)
plot(t, tau, 'r')
xlabel('Time')
ylabel('Control Torque')
title('Control Torque vs. Time')

% 绘制单摆小球的运动轨迹
figure(3)
x0 = 0; y0 = 0;
v = VideoWriter('pendulum.avi');
open(v);

for i = 1:length(t)
    theta = y(i, 1);
    x = l * sin(theta);
    y = -l * cos(theta);
    
    plot([x0 x], [y0 y], 'b', 'LineWidth', 2);
    hold on;
    plot(x, y, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
    axis([-1.5 1.5 -1.5 0.5])
    axis equal
    grid on
    
    frame = getframe(gcf);
    writeVideo(v, frame);
    clf;
end

close(v);

% 微分方程组
function dy = ode_pendulum(t, y)
    global m l g
    dy = zeros(2, 1);
    
    % PID控制器
    theta_desired = pi/6;
    theta_error = theta_desired - y(1);
    dy(1) = y(2);
    dy(2) = -3 * g / (m * l) * sin(y(1)) + 3 * (Kp * theta_error + Ki * trapz(theta_error) + Kd * y(2));
end

运行完上述代码后，会生成三个图像，分别是角度随时间变化的曲线、控制力矩随时间变化的曲线和单摆小球的运动轨迹。同时，也会生成一个名为“pendulum.avi”的视频文件，保存了单摆小球的模拟实物动态图。

需要注意的是，由于每个计算机的性能不同，可能需要对PID控制器的参数进行调整，才能够得到较好的控制效果。