编辑
2025-04-20
自动控制
00

目录

论文
数学建模公式
1. 参数定义
2. 动能(K)
3. 势能(V)
4. 拉格朗日方程
5. 方程推导
Brushed DC Motor Parameter
有刷直流电机参数
电气部分动态方程
机械部分动态方程
Controllers
Sound Generation
Compare PID controller and LQR controller
Conclusion
Acknowledgments
Reference

《SWING UP AND BALANCING OF A REACTION WHEEL INVERTED PENDULUM》

https://github.com/B-Paweekorn/Reaction-wheel-inverted-pendulum

论文

在定义了坐标之后,下一步是用广义坐标计算动能和总势能,然后,由于能量是标量函数,只需将动能和势能相减即可得到系统中的总机械能。通过以下表达式可以找到动能和势能的拉格朗日量LL

L=TVL = T - V

其中,TT是动能,VV是势能。

由于拉格朗日量是广义坐标及其导数的函数,运动方程如下所示:

ddt(Lq˙k)Lqk=τk\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_k}\right) - \frac{\partial L}{\partial q_k} = \tau_k

其中,τk\tau_{\boldsymbol{k}}是沿qkq_{\boldsymbol{k}}方向的广义力或力矩。

2.2 状态空间建模

状态空间是一种数学模型,用于将物理系统表示为输入、输出和状态的关系。状态空间模型的一般形式由以下方程给出:

X˙=AX+BUY=CX+DU\begin{aligned} & \dot{X} = A X + B U \\ & Y = C X + D U \end{aligned}

其中:
XX 是状态向量,

YY 是输出向量,

UU 是输入向量,

AA 是状态矩阵,

BB 是输入矩阵,

CC 是输出矩阵,

DD 是前馈矩阵,

X˙\dot{X} 是状态向量的一阶导数。

2.3 电机模型

直流(DC)电机的机械系统数学模型通常基于电机转矩进行描述。在本实验中,直流电机用于驱动反作用轮。在直流电机中,电机转矩与电机电流成正比。由于我们使用 PWM(脉宽调制)电压驱动器而非电流驱动器来驱动电机,因此需要电机的动态模型来控制和近似转矩。

image.png

通常,直流电机产生的转矩(τ\boldsymbol{\tau})与转子线圈中的电流成正比,可以用直流电机的机电方程表示:

τ=kti其中,kt 是电机的转矩常数。\tau = k_t \cdot i \quad \text{其中,} \quad k_t \text{ 是电机的转矩常数。}

反电动势(e\boldsymbol{e})与电机转子的角速度成正比,可以用直流电机的电气方程表示:

e=keθ˙其中,ke 是电机的反电动势常数,θ˙ 是电机轴的角速度。\begin{aligned} & \boldsymbol{e} = k_e \cdot \dot{\theta} \quad \text{其中,} \quad k_e \text{ 是电机的反电动势常数,} \\ & \dot{\theta} \text{ 是电机轴的角速度。} \end{aligned}

假设直流电机的效率为 100%100\%,则转矩常数 ktk_t 和反电动势常数 kek_e 应相等。根据基尔霍夫电压定律,电机电路可以推导出以下表达式:

vRiLdidtkeθ˙=0v - R i - L \frac{d i}{d t} - k_e \dot{\theta} = 0

2.4 起摆控制器(Bang-Bang控制器)

该系统的起摆控制采用Bang-Bang控制器,其输出信号在输入信号的上限和下限之间切换。Bang-Bang控制器可以控制摆杆的运动,但无法精确控制摆杆在接近直立位置时的速度,而这一点对稳定控制至关重要。因为如果摆杆到达直立位置时的角速度过高,反作用轮的扭矩可能不足以使其停止。因此,当摆杆接近直立位置时,需要切换至另一种控制器(如LQR控制器)以实现稳定控制。

Bang-Bang控制器的激活由一个基于能量的控制策略决定。在Bang-Bang控制器增加摆杆能量的同时,能量控制器会实时检测摆杆的机械能量,并将其与目标能量进行比较。如果当前能量未达到目标值,则Bang-Bang控制器保持开启;否则,Bang-Bang控制器关闭,电机停止工作。

image.png

2.5 线性二次调节器(LQR控制器)

LQR是一种用于优化线性时不变系统控制的数学框架,其目标是通过最小化一个成本函数(通常表示为控制输入随时间的变化)来实现最优控制。该控制器基于当前系统状态计算反馈控制律,以确保系统稳定性并达成控制目标。

2.6 PID控制器

PID是“比例-积分-微分控制器”的缩写,是动态控制系统中应用最广泛的控制器。其控制输出由比例项、积分项和微分项叠加而成,表达式如下:

u(t)=Kpe(t)+Ki0τe(τ)dτ+Kdddte(t)u(t) = K_p e(t) + K_i \int_0^\tau e(\tau) d\tau + K_d \frac{d}{dt} e(t)

其中:
u(t)u(t) 为控制输入信号,

KpK_p 为比例增益,

KiK_i 为积分增益,

KdK_d 为微分增益,

e(t)e(t) 为误差信号,

tt 为瞬时时间,

τ\tau 为积分变量。

本实验将结合LQR控制器与PID控制器共同调节摆杆角速度,以对比两种控制器的输出性能与稳定性差异。

  1. 方法论

在本项目框架中,系统框图由两大核心部分组成:控制器部分(包含关键的起摆控制器和平衡控制器)以及作为重要测量工具的仿真部分。仿真部分进一步划分为两个模块:反应轮倒立摆动力学模型和反应轮倒立摆仿真模型。这些组件共同构成了一个用于控制策略分析与验证的完整框架。

image.png

数学建模公式

1. 参数定义

L1: 摆杆质心到转轴的长度。

L2: 飞轮安装位置到转轴的长度。

m1, m2: 摆杆和飞轮的质量。

θp, θr: 摆杆和飞轮的角位移(飞轮角为相对摆杆的相对角)。

I1, I2: 摆杆和飞轮(含电机)的转动惯量。

Tm: 电机施加的扭矩(驱动飞轮)。

Td: 外部扰动扭矩。

2. 动能(K)

系统动能包含两部分:

  1. 摆杆的平动与转动动能:
    K1=12(m1L12+m2L22+I1)θ˙p2K_1 = \frac{1}{2}(m_1 L_1^2 + m_2 L_2^2 + I_1)\dot{\theta}_p^2
  2. 飞轮的转动动能:
    飞轮的总角速度为 θ˙p+θ˙r\dot{\theta}_p + \dot{\theta}_r,因此:
    K2=12I2(θ˙p+θ˙r)2K_2 = \frac{1}{2}I_2(\dot{\theta}_p + \dot{\theta}_r)^2

展开后总动能为:

K=12(m1L12+m2L22+I1+I2)θ˙p2+I2θ˙pθ˙r+12I2θ˙r2K = \frac{1}{2}(m_1 L_1^2 + m_2 L_2^2 + I_1 + I_2)\dot{\theta}_p^2 + I_2 \dot{\theta}_p \dot{\theta}_r + \frac{1}{2}I_2 \dot{\theta}_r^2

3. 势能(V)

摆杆和飞轮的势能由其质心高度决定:

V=(m1L1+m2L2)gcosθpV = (m_1 L_1 + m_2 L_2)g \cos\theta_p

4. 拉格朗日方程

拉格朗日量 L=KVL = K - V,对广义坐标 θr\theta_rθp\theta_p 分别应用方程:

ddt(Lθ˙r)Lθr=Qr,ddt(Lθ˙p)Lθp=Qp\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}_r}\right) - \frac{\partial L}{\partial \theta_r} = Q_r, \quad \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}_p}\right) - \frac{\partial L}{\partial \theta_p} = Q_p

其中广义力: • Qr=TmQ_r = T_m(电机扭矩作用于飞轮)。

Qp=Tm+TdQ_p = -T_m + T_d(摆杆受到的反作用扭矩和扰动)。

L=KV=[12(m1L12+m2L22+I1+I2)θ˙p2+I2θ˙pθ˙r+12I2θ˙r2](m1L1+m2L2)gcosθpL = K - V = \left[ \frac{1}{2}(m_1 L_1^2 + m_2 L_2^2 + I_1 + I_2)\dot{\theta}_p^2 + I_2 \dot{\theta}_p \dot{\theta}_r + \frac{1}{2}I_2 \dot{\theta}_r^2 \right] - (m_1 L_1 + m_2 L_2)g \cos\theta_p

5. 方程推导

θp\theta_p 的方程

L=KV=[12(m1L12+m2L22+I1+I2)θ˙p2+I2θ˙pθ˙r+12I2θ˙r2](m1L1+m2L2)gcosθpL = K - V = \left[ \frac{1}{2}(m_1 L_1^2 + m_2 L_2^2 + I_1 + I_2)\dot{\theta}_p^2 + I_2 \dot{\theta}_p \dot{\theta}_r + \frac{1}{2}I_2 \dot{\theta}_r^2 \right] - (m_1 L_1 + m_2 L_2)g \cos\theta_p
  1. 计算导数项:
    Lθ˙p=(m1L12+m2L22+I1+I2)θ˙p+I2θ˙r\frac{\partial L}{\partial \dot{\theta}_p} = (m_1 L_1^2 + m_2 L_2^2 + I_1 + I_2)\dot{\theta}_p + I_2 \dot{\theta}_r
    ddt(Lθ˙p)=(m1L12+m2L22+I1+I2)θ¨p+I2θ¨r\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}_p}\right) = (m_1 L_1^2 + m_2 L_2^2 + I_1 + I_2)\ddot{\theta}_p + I_2 \ddot{\theta}_r
  2. 势能导数项: ( cos ⁡ x ) ′ = − sin ⁡ x
    Lθp=(m1L1+m2L2)gsinθp\frac{\partial L}{\partial \theta_p} = (m_1 L_1 + m_2 L_2)g \sin\theta_p
  3. 代入拉格朗日方程:
    (m1L12+m2L22+I1+I2)θ¨p+I2θ¨r(m1L1+m2L2)gsinθp=Tm+Td(2)(m_1 L_1^2 + m_2 L_2^2 + I_1 + I_2)\ddot{\theta}_p + I_2 \ddot{\theta}_r - (m_1 L_1 + m_2 L_2)g \sin\theta_p = -T_m + T_d \tag{2}

θr\theta_r 的方程

L=KV=[12(m1L12+m2L22+I1+I2)θ˙p2+I2θ˙pθ˙r+12I2θ˙r2](m1L1+m2L2)gcosθpL = K - V = \left[ \frac{1}{2}(m_1 L_1^2 + m_2 L_2^2 + I_1 + I_2)\dot{\theta}_p^2 + I_2 \dot{\theta}_p \dot{\theta}_r + \frac{1}{2}I_2 \dot{\theta}_r^2 \right] - (m_1 L_1 + m_2 L_2)g \cos\theta_p
  1. 计算导数项:
Lθ˙r=Kθ˙r=I2(θ˙p+θ˙r)\frac{\partial L}{\partial \dot{\theta}_r} = \frac{\partial K}{\partial \dot{\theta}_r} = I_2 (\dot{\theta}_p + \dot{\theta}_r)
ddt(Lθ˙r)=I2(θ¨p+θ¨r)\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\theta}_r} \right) = I_2 (\ddot{\theta}_p + \ddot{\theta}_r)
Lθr=0\frac{\partial L}{\partial \theta_r} = 0
  1. 代入拉格朗日方程:
    I2(θ¨p+θ¨r)=Tm(1)I_2 (\ddot{\theta}_p + \ddot{\theta}_r) = T_m \tag{1}

DC Motor Dynamics

image.png

Brushed DC Motor Parameter

有刷直流电机参数

Vin - 输入电压

R - 电机电阻

L - 电机电感

i - 电机电流

B - 电机阻尼系数

J - 电机转子转动惯量

ke - 反电动势常数

kt - 扭矩常数

电气部分动态方程

假设电感效应可忽略(LRL \ll R),电压平衡方程简化为:

Vin=Ri+keθ˙rV_{\text{in}} = R \, i + k_e \dot{\theta}_r

由此可得电枢电流 ii 的表达式:

i=Vinkeθ˙rRi = \frac{V_{\text{in}} - k_e \dot{\theta}_r}{R}

其中:
θ˙r\dot{\theta}_r 为电机转子角速度

keθ˙rk_e \dot{\theta}_r 表示反电动势电压

机械部分动态方程

电机扭矩 TmT_m 与电枢电流 ii 的关系由扭矩常数 ktk_t 决定:

Tm=kti=ktVinkeθ˙rRT_m = k_t \, i = k_t \frac{V_{\text{in}} - k_e \dot{\theta}_r}{R}

扭矩平衡方程为:

Tm=Jθ¨r+Bθ˙rT_m = J \ddot{\theta}_r + B \dot{\theta}_r

其中:
θ¨r\ddot{\theta}_r 为电机转子角加速度

Bθ˙rB \dot{\theta}_r 表示阻尼扭矩

该推导表明,电机扭矩 TmT_m 同时受电气输入和机械运动的影响。

Controllers

  • LQR Controller: A linear quadratic regulator designed to stabilize the pendulum in the upright position.

    Linearization dynamics model

    • Part RWIP

      when sinθp -> 0 sinθp = θp

        θp¨=m1gL1θp + m2gL2θp  Tm+ Tdm1L12+m2L2+I1\ddot{\theta_p} = \frac{m_{1}gL_{1}\theta_{p}\ +\ m_{2}gL_{2}\theta_{p}\ -\ T_{m} +\ T_{d}}{m_{1}L_{1}^{2}+m_{2}L_{2}+I_{1}}

        θr¨=TmI2\-m1gL1 + m2gL2  Tm+ Tdm1L12+m2L2+I1\ddot{\theta_r} = \frac{T_{m}}{I_{2}}\-\frac{m_{1}gL_{1}\ +\ m_{2}gL_{2}\ -\ T_{m} +\ T_{d}}{m_{1}L_{1}^{2}+m_{2}L_{2}+I_{1}}

    • Part Motor

       We can estimate that L << R

        Vin=Ri+keθrVin = R i + k_e θ_r

        Tm=ktiT_{m} = k_t i

    State space

    The proceeding equations are valid around the operating point where θp = 0

[θ˙pθ¨pθ˙rθ¨r]=[0100(m1L1+m2L2)gm1L12+m2L22+J00ktke(m1L12+m2L22+J)R0001(m1L1+m2L2)gm1L12+m2L22+J00(m1L12+m2L22+2J(m1L12+m2L22+J)J)(ktkcR)][0θpθ˙pθrθ˙r]+[ktm1L12+m2L22+J)R0(m1L12+m2L22+2J(m1L12+m2L22+J)J)(ktR)]Vin\left[\begin{array}{c} \dot{\theta}_p \\ \ddot{\theta}_p \\ \dot{\theta}_r \\ \ddot{\theta}_r \end{array}\right]=\left[\begin{array}{cccc} 0 & 1 & 0 & 0 \\ \frac{\left(m_1 L_1+m_2 L_2\right) g}{m_1 L_1^2+m_2 L_2^2+J} & 0 & 0 & \frac{k_t k_e}{\left(m_1 L_1^2+m_2 L_2^2+J\right) R} \\ 0 & 0 & 0 & 1 \\ -\frac{\left(m_1 L_1+m_2 L_2\right) g}{m_1 L_1^2+m_2 L_2^2+J} & 0 & 0 & -\left(\frac{m_1 L_1^2+m_2 L_2^2+2 J}{\left(m_1 L_1^2+m_2 L_2^2+J\right) J}\right)\left(\frac{k_t k_c}{R}\right) \end{array}\right]\left[\begin{array}{c} 0 \\ \theta_p \\ \dot{\theta}_p \\ \theta_r \\ \dot{\theta}_r \end{array}\right]+\left[\begin{array}{c} k_t \\ \frac{\left.m_1 L_1^2+m_2 L_2^2+J\right) R}{0} \\ \left(\frac{m_1 L_1^2+m_2 L_2^2+2 J}{\left(m_1 L_1^2+m_2 L_2^2+J\right) J}\right)\left(\frac{k t}{R}\right) \end{array}\right] V_{i n}
  • PID Controller: The stabilize controller to compare with LQR

    Transfer function

        θp(s)τm(s)=sJm2L22s3+(BI1+B+dpJ+m2L22)s2((m1L1+m2L2)g(J+m2L22)I1Bdp(J+m2L22)I1)s(m1L1+m2L2)Bg(J+m2L22)I1\Large\frac{\theta_{p}(s)}{\tau_{m}(s)}=\frac{\frac{s}{-J-m_{2}L_{2}^{2}}}{s^{3}+\left(\frac{B}{I_{1}}+\frac{B+d_{p}}{J+m_{2}L_{2}^{2}}\right)s^{2}-\left(\frac{\left(m_{1}L_{1}+m_{2}L_{2}\right)g}{\left(J+m_{2}L_{2}^{2}\right)I_{1}}-\frac{B\cdot d_{p}}{\left(J+m_{2}L_{2}^{2}\right)I_{1}}\right)s-\frac{\left(m_{1}L_{1}+m_{2}L_{2}\right)Bg}{\left(J+m_{2}L_{2}^{2}\right)I_{1}}}

    Root Locus Design

    To predict the system's characteristics as the gain (Kp) is adjusted and poles move, design the root locus.

        image

    Root Locus of the system by default parameter set. It has one zero (s = 0), and three poles (s = -74.93; s =-3.88e-3; s = 73.53).

    Closed Loop Root Locus

        image

      Where G(s) is Kp the gain can be adjusted ti make the closed loop poles to be in stable location The resultant Root Locus can be seen below (note to plot this graph in param.py you need to set Stabilize_Controller to "PID" mode and set plot_rootlocus to "True") in this graph you can click pole position you want to know Gain Kp to adjust your system characteristics.

       image

  • Bang-bang Controller: The swing up control routine and the stabilizing control routine are switched between -25 to 25 degree

  

  • Brake Controller: Used as reduced energy of RWIP when RWIP have too much energy for stabilze

Sound Generation

        The simulation incorporates sound generation related to the speed of the reaction wheel. This feature adds an auditory element to the simulation, enhancing the user experience.


Compare PID controller and LQR controller

In this project, we explore and compare the performance of two different control strategies: PID (Proportional-Integral-Derivative) controller and LQR (Linear Quadratic Regulator) controller.
image

Linear quadratic regulator

CodeCogsEqn

Error (deg)settling time (s)Power (Watt)
50.660.6
60.731.01
70.851.75
81.073.5
9can't stabilizecan't stabilize

Max Disturbance : 9.32 Nm

PID : Kp = 500 (choose form root locus)

Error (deg)settling time (s)Power (Watt)
520.4310.06
621.0212.17
721.3713.95
821.5915.96
9can't stabilizecan't stabilize

Max Disturbance : 8.05 Nm

PID : Kp = 215800 (choose form root locus) image Notice that when choose unstable pole the system still stable because now it have hardware limit so the character of controller same like Fuzzy logic control to see unstable you need to unlock hardware limitation by set MotorLimit to False in param.py

Conclusion

Stabilize boundary
LQRCan stabilize in every position
PIDCan stabilize only in small boundary

PID Controller The PID controller is a widely used feedback control system that relies on three components: Proportional, Integral, and Derivative. Here's a brief overview of each component:

  • Proportional (P): Reacts to the current error.
  • Integral (I): Reacts to the accumulation of past errors.
  • Derivative (D): Predicts future errors based on the rate of change.

Advantages of PID:

  • Simplicity and ease of implementation.
  • Effectiveness in a wide range of systems.
  • Intuitive tuning parameters for performance optimization.

Considerations:

  • Tuning may be required for optimal performance in different systems.
  • Limited capability to handle complex or nonlinear systems.

LQR Controller The LQR controller is designed based on the principles of optimal control theory. It minimizes a cost function that combines both state and control input, making it suitable for linear, time-invariant systems.

Advantages of LQR:

  • Optimal control solution for linear systems.
  • Ability to handle systems with multiple inputs and outputs.
  • Incorporates a mathematical model for optimal performance.

Considerations:

  • Requires a good understanding of the system dynamics for effective modeling.
  • Limited applicability to strictly linear systems.

Acknowledgments

        This project is part of the coursework for FRA333 Robot Kinematics at the Institute of Field Robotics, King Mongkut’s University of Technology Thonburi. Special thanks to the course instructors for their guidance and support.

Feel free to explore, modify, and extend this project for educational and research purposes.


Reference

如果对你有用的话,可以打赏哦
打赏
ali pay
wechat pay

本文作者:Dong

本文链接:

版权声明:本博客所有文章除特别声明外,均采用 CC BY-NC。本作品采用《知识共享署名-非商业性使用 4.0 国际许可协议》进行许可。您可以在非商业用途下自由转载和修改,但必须注明出处并提供原作者链接。 许可协议。转载请注明出处!