整体撸一遍PMSM的滑模观测器(SMO)公式+模型+代码

文章介绍了永磁同步电机无速度传感器控制中滑模观测器的设计原理和稳定性分析。通过电流模型建立观测器,并利用李亚普诺夫稳定性理论证明了观测器的稳定性。在仿真模型的搭建和参数调试中,探讨了滑模增益和低通滤波器对观测效果的影响,并解决了角度提取过程中的问题,实现了转速闭环控制。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

 滑模的基本原理参照这一系列文章:

 滑模系列文章链接:

永磁同步电机矢量控制到无速度传感器控制学习教程(PMSM)(一)

永磁同步电机矢量控制基础补充(五)——什么是低通滤波器?

永磁同步电机无速度传感器控制(一)——滑模观测器(一)【位置估计原理】

永磁同步电机无速度传感器控制(一)——滑模观测器(二)【滑模观测器设计过程】

永磁同步电机无速度传感器控制(一)——滑模观测器(三)【由扩展反电势得到电机位置和速度信息】

永磁同步电机无速度传感器控制(一)——滑模观测器(四)【仿真搭建及其结果分析】

该文章中包含其他滑模观测器更为详细部分的讲解链接。在进入代码学习之前,可以先对这个部分进行阅读与学习。

本篇文章意图记录模观测器代码学习和调试过程。在此先仅考虑表贴式(隐极)电机,即Ld = Lq。

1、电机公式模型及观测器的设计原理

根据前面的一些基础的理解,滑模观测器观测器最终得到角度的方法是对反电动势 Ealpha 和 Ebeta 求取反正切得到的,因此需要分析的数学模型将从矢量控制中dq坐标系转到 alpha-beta 两相静止坐标系中,两相静止坐标系的电压方程如下所示:

 其中反电动势的表达式如下所示:(注意这些分析都是基于表贴式永磁同步电机Ld=Lq,需要研究凸极电机的这个公式并不适用)

 由于电压一般进作为激励物理量,而电流才是响应物理量,另外由于电机电流模型中包含常微分方程形式(y'  = f (x,y) ),这种形式的观测器结构更利于稳定性的推导与证明。因此电机中观测器的设计一般是基于电流模型。博客中介绍的MARS其实也是一种电流模型,其稳定性推导就是基于常微分方程展开的popov稳定性理论,滑模也选择这个模型进行。为此将上面两个公式所示的电压模型进行一定的公式变化,得到如下图所示的电流方程:

 依据电机模型,等阶次的、格式对应的建立滑模观测器模型,如下所示:

 其中观测所得的反电动势表达式如下图所示:

误差方程表达式为(用估测模型 - 实际电机模型)

2、滑模观测器的稳定性分析

一个观测器能否渐进稳定是决定观测器是否能够使用的关键,上面虽然设计了滑模观测器结构,但是它具体能不能够正常使用,还是一个未知的问题。因此有必要对齐进行稳定性证明。稳定性的证明需要用到李雅普诺夫稳定性理论中的稳定性判据。

李亚普诺夫稳定性理论基于李亚普诺夫函数的概念,该函数是一个实值函数,能够度量系统状态的变化情况。根据李亚普诺夫函数的性质,可以推导出系统的稳定性条件。如果李亚普诺夫函数对于系统的所有状态都是非正的,并且在系统中存在一个李亚普诺夫函数值为零的稳定点,那么系统就被认为是稳定的。

简单来说,如果平衡状态 x 受到扰动后,仍然停留在 x 附近,我们就称 x 在李雅普诺夫意义下是稳定的(Lyapunov stable)。

更进一步,如果平衡状态 xe 受到扰动后,最终都会收敛到 x ,我们就称 x 在李雅普诺夫意义下是渐进稳定的(Asymptotically stable)。

再进一步,如果平衡状态 x 受到任何扰动后,最终都会收敛到 x ,我们就称 x 在李雅普诺夫意义下是大范围内渐进稳定的(Asymptotically stable in large)

由于滑模是基于开关函数的,因此滑模观测器并不属于一种线性系统,属于非线性系统。因此应用李亚普诺夫第二法较为合适(李亚普诺夫第一方法需要讨论线性化系统的特征值分布来研究原非线性系统的稳定性问题,这是一种间接的证明方法也称间接法)。李亚普诺夫第二法的判据为:

  1. 李亚普诺夫第二法则:对于连续时间系统,如果存在一个连续可微的李亚普诺夫函数V(x),满足以下两个条件:

    • V(0) = 0,当且仅当 x = 0 时成立。
    • V(x) > 0,对于所有 x ≠ 0 成立。
    • 对于系统的状态变化,V(x) 的导数满足 dV(x)/dt < 0,即 V(x) 随时间递减。

    则系统是渐近稳定的。总结成方程形式就是  

通过李亚普诺夫第二法则,能够证明稳定的系统能量总是不断被耗散的,李雅普诺夫通过定义一个标量函数 V(x)(通常能代表广义能量)来分析稳定性。这种方法避免了复杂的方程求解和线性化,直接证明系统的稳定性,因此称为直接法。对于滑模观测器(其他观测器也类似),稳定判据中的V(x)其实就是误差函数。因为误差函数代表了系统的熵增和熵减,对应李亚普诺夫函数的能量耗散原理。

误差函数为:

 根据电流方程及电流误差方程:

可以推出:

β轴的与此同理,满足下面条件即可证明稳定性:

由于平方项必定大于0,因此得到如下推导:

 最终得到系统稳定的条件为

 这样就能产生滑动运动,这样误差动态方程就是渐进稳定的,也就保证了滑模观测器的全局稳定性,用起来就不担心。

3、连续域仿真模型的搭建

3.1 滑模前期的仿真基础

依据上述阐述的原理对仿真模型进行搭建,验证算法的有效性。首先需要准备一个基础的双闭环仿真,先尝试用实际速度做闭环时,滑模是否能够估计到真实转速。电机参数如下图所示。

3.2 搭建滑模观测器 

在此基础上搭建滑模观测器,避免大家上面再去翻,再次粘贴一次公式:

 3.3 滑模参数调试

滑模观测器中除了电机参数需要填入外,还需要填入两个关键参数——滑模增益和低通滤波器截止频率。由于本文所使用的电机母线电压只有24V,反电动势比不可能大于24,根据前文稳定性推导部分可知,设置滑模增益为24系统必定处于全局可收敛。本仿真系统的开关频率时10kHz,电机额定频率为250Hz,因此考虑设置滑模后端的低通滤波器为中间值即1000Hz(1000Hz对应的wc = 1000*2*pi = 6280)。带入上述参数进行调试,下图为观测所得Ealpha和Ebeta波形。

 根据上图所示波形,可知该滑模观测器已基本实现反电动势波形的观测,但是明显可以看到此时的波形毛刺较大,抖振较为严重,考虑降低滑模增益调试。初始给定的24V仅仅是一个较为保守的值,将其缩小根号3倍,即 24/1.732 = 13.85。该滑模增益时反电动势波形如下图所示,此时抖振已被减小的较小了。

 此时仍然存在的一些谐波是滑模观测器自身在上下滑动过程中造成的,并非滑模增益设置不够合适,而是低通滤波器的截止频率可能太高了,许多谐波没有被滤掉。因此考虑降低低通滤波器的截止频率到2倍额定频率 = 250Hz*2 = 500Hz,将其转化为wc = 500*2*pi = 3140。反电动势波形如下图所示,其实调到如下程度就可以了,可进行转速和位置的提取。

3.3 袁雷书上SMO_acrtan模型存在的一点遗漏

        刚开始看书的时候发现这个地方怎么也看不懂,经过大佬的指导发现这里存在一些问题,(感谢千羽同学提供的一些思路)将其总结如下面4个,后续3.4部分探究解决过程。具体的如下图所示,袁雷书上模型的链接如下:

链接:https://pan.baidu.com/s/1FIrezc-PAqBylqqVvO0IaA 
提取码:8888

 3.4 重新推了一下基于反正切函数的转子位置角度提取

根据电机反电动势表达式,对反电动势求反正切即可得到电机的角度信息。

 电机转子位置表达如下所示:

 搭建仿真实现其过程,仿真模块搭建如下图所示:

         此时将电机转子位置角与滑模估算角度放在一个示波器中进行观察,如下图所示,此时波形非常的奇怪,电机位置角(红色)一直在累加,这是因为simulink仿真自带的PMSM模块出来的电角度就是会一直叠加,为此我们需要在其前级加上一个取余函数。编写一个s-function或者搭建一个取余模块,为什么这里写出来呢,因为我在这里被坑了好久,有的版本找不到mod模块,最好用s-function写一个mod()函数的处理一下。

         最终得到的转子位置对比图如下图所示:

         由于arctan函数的输出值(-pi/2,pi/2),其函数图像如下图所示。

        需要具体的分析 arctan函数输出值与实际电机位置角的关系,当acrtan(-Ealpha /Ebeta) = 1时,实际存在两个可能的结果,一个是pi/4,一个是5pi/4,因为sin(pi/4) = 0.707 ,cos(pi/4)=0.707,但是sin(5pi/4)=-0.707,cos(5pi/4)=-0.707。因此仅凭acrtan(-Ealpha /Ebeta) = 1无法确定具体是pi/4还是5pi/4,根据一系列推论总结可得(下面这个公式其实就是通过每个象限的值一一枚举然后总结出来的,就把1 2 3 4 象限的不同的值和 arctan 的输出值相互对应出来的)

其模块代码为:

function y = fcn(Ealpha,Ebeta,Acrtan)
%#codegen
x = 0;
if Acrtan > 0 
    if Ebeta < 0
        x = pi + Acrtan;
    elseif  Ebeta > 0
        x =  Acrtan;
    elseif Ebeta == 0
        x = 0.5*pi;
    end
elseif Acrtan<0
    if Ebeta < 0
        x = pi + Acrtan;
    elseif  Ebeta > 0
        x = 2*pi + Acrtan;
    elseif Ebeta == 0
        x = 1.5*pi;
    end
elseif Acrtan == 0
     if Ebeta < 0
        x = pi;
    elseif  Ebeta > 0
        x = 0;
     elseif Ebeta == 0
        x = 0;
    end       
end
y = x;

 最终结果对比,蓝色为滑模观测器结果,黄色为电机实际角度,可以看到基本一致,由于低通滤波器的存在,有一定的延迟。

     低通滤波器造成的延迟会影响电机转子位置的准确性,因此需要对这个延迟进行一定的补偿,其实低通滤波器的截止频率我们是知道的,补偿这个延迟理论上是较好进行的。j将低通滤波器的传递函数转化为幅频和相频公式可得:

         将滑模所用的低通滤波器传递函数转换一下可得:

         其相位延迟等于:

         最终得到补偿相位公式为:

         将其补偿到输出电机角度上后结果如下图所示,估算角度基本与实测电机角度重合。但是此时我的补偿角度是拿电机实际we算出来的,所以波形效果相当好。(注意这里如果直接补偿上去,会导致出现大于2pi的角度,此时就需要稍微处理一下,把超出2pi的角度转为小角度,也就是减去2pi即可)

         实际是需要通过反电动势去算取we,其计算公式为:

         最终观测输出结果如下图所示,效果也是挺好的。 

整体获取角度+补偿过程的仿真截图如下图所示

最终实现转速闭环控制效果如下图所示:

四、滑模观测器的离散化

        仿真中可以按照连续模型去搭建,在代码中则必须考虑MCU的离散环境,先看连续域的电机电流方程,如下图所示。式子中 R 代表定子电阻,L 为定子电感,Valpha 为两相静止坐标系下的alpha 电压,ealpha 为两相静止坐标系下的 alpha 反电动势。根据 ealpha 和 ebeta 的表达式可知假如能够观测反电动势的大小,对齐求反正切即可得到电机电角度信息。

依据前向欧拉方程进行离散化推导,离散化操作可参照这篇文章。

数字滤波器的实现——低通滤波器再探究_低通滤波器差分方程_沉沙丶的博客-CSDN博客

        

以下是一个基于STM32的PMSM电机控制器中,使用了基于滑模观察器(SMO)和锁相环(PLL)的代码示例: ``` #include "stm32f4xx.h" #include "arm_math.h" #define PI 3.14159265358979323846f #define SQRT3 1.73205080756887729352f // Motor parameters (需要根据具体电机的参数进行修改) #define R_PHASE 0.83f #define L_D 0.000354f #define L_Q 0.000354f #define POLE_PAIRS 7 // SMO parameters #define SMO_LAMBDA 10.0f #define SMO_BETA 1.0f // PLL parameters #define PLL_KP 0.5f #define PLL_KI 0.01f // Sampling frequency #define FS 10000.0f // Global variables float32_t V_alpha, V_beta, V_d, V_q, I_alpha, I_beta, I_d, I_q, theta_elec, theta_mech; float32_t Id_error, Iq_error, Id_ref, Iq_ref, Id_out, Iq_out, Vd_out, Vq_out; float32_t smo_estimated_flux_d, smo_estimated_flux_q, smo_estimated_speed; float32_t pll_estimated_speed, pll_error, pll_integral, pll_proportional; float32_t sector_angle, sector_duty_cycle, period_ticks, pwm_ticks; float32_t sin_theta, cos_theta, sin_theta_120, cos_theta_120, sin_theta_240, cos_theta_240; // Sine and cosine lookup tables const float32_t sin_table[360] = {...}; const float32_t cos_table[360] = {...}; // SMO function void SMO(float32_t V_alpha, float32_t V_beta, float32_t I_alpha, float32_t I_beta, float32_t dt) { float32_t Lm, Ls, Rs, inv_Ls, sm_alpha, sm_beta, sm_d, sm_q, sm_flux_norm; float32_t smo_error_d, smo_error_q, smo_error_norm, smo_psi_alpha, smo_psi_beta; Lm = L_D + L_Q + (R_PHASE * dt) / 2.0f; Ls = L_D + L_Q - (R_PHASE * dt) / 2.0f; Rs = R_PHASE; inv_Ls = 1.0f / Ls; // Clarke transform sm_alpha = I_alpha; sm_beta = SQRT3 * I_beta - SQRT3 / 2.0f * I_alpha; // Park transform sm_d = cos_theta * sm_alpha + sin_theta * sm_beta; sm_q = -sin_theta * sm_alpha + cos_theta * sm_beta; // Compute estimated flux smo_estimated_flux_d += dt * (-smo_estimated_flux_d * SMO_LAMBDA / Lm + sm_d); smo_estimated_flux_q += dt * (-smo_estimated_flux_q * SMO_LAMBDA / Lm + sm_q); sm_flux_norm = arm_sqrt_f32(smo_estimated_flux_d * smo_estimated_flux_d + smo_estimated_flux_q * smo_estimated_flux_q); // Compute error between estimated and actual flux smo_psi_alpha = Lm * Id_out + smo_estimated_flux_d; smo_psi_beta = Lm * Iq_out + smo_estimated_flux_q; smo_error_d = sm_d - (smo_psi_alpha * cos_theta + smo_psi_beta * sin_theta) * inv_Ls; smo_error_q = sm_q - (-smo_psi_alpha * sin_theta + smo_psi_beta * cos_theta) * inv_Ls; smo_error_norm = arm_sqrt_f32(smo_error_d * smo_error_d + smo_error_q * smo_error_q); // Update estimated speed smo_estimated_speed = (smo_estimated_flux_q * smo_error_d - smo_estimated_flux_d * smo_error_q) / (sm_flux_norm * sm_flux_norm + SMO_BETA); // Compute d and q current estimates Id_out = smo_psi_alpha * inv_Ls + smo_estimated_speed * smo_estimated_flux_q / sm_flux_norm; Iq_out = smo_psi_beta * inv_Ls - smo_estimated_speed * smo_estimated_flux_d / sm_flux_norm; } // PLL function void PLL(float32_t V_alpha, float32_t V_beta, float32_t I_alpha, float32_t I_beta, float32_t dt) { float32_t V_d_filtered, V_q_filtered, I_d_filtered, I_q_filtered, V_cross, I_cross; // Clarke transform I_alpha = I_alpha; I_beta = SQRT3 * I_beta - SQRT3 / 2.0f * I_alpha; // Park transform I_d = cos_theta * I_alpha + sin_theta * I_beta; I_q = -sin_theta * I_alpha + cos_theta * I_beta; // Compute Vd and Vq V_d = V_alpha - Rs * I_d - smo_estimated_speed * smo_estimated_flux_q; V_q = V_beta - Rs * I_q + smo_estimated_speed * smo_estimated_flux_d; // Low-pass filter Vd and Vq V_d_filtered += dt * (V_d - V_d_filtered) / (L_D * 2.0f * PI * 1000.0f + dt); V_q_filtered += dt * (V_q - V_q_filtered) / (L_Q * 2.0f * PI * 1000.0f + dt); // Low-pass filter Id and Iq I_d_filtered += dt * (I_d - I_d_filtered) / (L_D * 2.0f * PI * 1000.0f + dt); I_q_filtered += dt * (I_q - I_q_filtered) / (L_Q * 2.0f * PI * 1000.0f + dt); // Compute V and I cross products V_cross = V_d_filtered * I_q_filtered - V_q_filtered * I_d_filtered; I_cross = I_d_filtered * smo_estimated_flux_q - I_q_filtered * smo_estimated_flux_d; // Compute estimated speed pll_error = V_cross / (V_cross * V_cross + I_cross * I_cross); pll_integral += pll_error * dt; pll_proportional = pll_error * PLL_KP; pll_estimated_speed = pll_proportional + pll_integral * PLL_KI; } // PWM function void PWM(float32_t duty_cycle) { if (duty_cycle < 0.0f) { duty_cycle = 0.0f; } else if (duty_cycle > 1.0f) { duty_cycle = 1.0f; } if (duty_cycle > 0.0f) { TIM1->CCR1 = (uint16_t)(duty_cycle * pwm_ticks); TIM1->CCR2 = (uint16_t)(duty_cycle * pwm_ticks); TIM1->CCR3 = (uint16_t)(duty_cycle * pwm_ticks); } else { TIM1->CCR1 = 0; TIM1->CCR2 = 0; TIM1->CCR3 = 0; } } // Main function int main(void) { // Initialize GPIO pins, timers, ADC, etc. // Main loop while (1) { // Read ADC values and convert to currents and voltages V_alpha = ...; V_beta = ...; I_alpha = ...; I_beta = ...; // Compute electrical and mechanical angles theta_elec += 2.0f * PI * (pll_estimated_speed / (POLE_PAIRS * FS)); if (theta_elec >= 2.0f * PI) { theta_elec -= 2.0f * PI; } theta_mech = theta_elec / POLE_PAIRS; // Compute sector angle and duty cycle sector_angle = theta_mech / (2.0f * PI / 3.0f); if (sector_angle < 0.0f) { sector_angle += 3.0f; } if (sector_angle >= 3.0f) { sector_angle -= 3.0f; } if (sector_angle < 1.0f) { sector_duty_cycle = sin_theta / 2.0f + cos_theta / SQRT3 * sin_theta_120 / 2.0f; } else if (sector_angle < 2.0f) { sector_duty_cycle = cos_theta / SQRT3 * sin_theta_120 / 2.0f - sin_theta / 2.0f; } else { sector_duty_cycle = -cos_theta / SQRT3 * sin_theta_120 / 2.0f - cos_theta / 2.0f; } // Compute period and PWM ticks period_ticks = (uint32_t)(SystemCoreClock / (POLE_PAIRS * FS)); pwm_ticks = (uint32_t)(period_ticks * sector_duty_cycle); // Run SMO and PLL SMO(V_alpha, V_beta, I_alpha, I_beta, 1.0f / FS); PLL(V_alpha, V_beta, I_alpha, I_beta, 1.0f / FS); // Compute Id and Iq references Id_error = Id_ref - Id_out; Iq_error = Iq_ref - Iq_out; // Run PI controllers for Id and Iq Vd_out = Id_error * L_D * 2.0f * PI * 1000.0f + Iq_error * smo_estimated_speed * L_Q * 2.0f * PI * 1000.0f; Vq_out = Iq_error * L_Q * 2.0f * PI * 1000.0f - Id_error * smo_estimated_speed * L_D * 2.0f * PI * 1000.0f; // Run PWM PWM(sector_duty_cycle); } } ```
评论 55
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值