【异步电机系列】电机参数离线辨识(含源码实现)

【一、闲话

        很久没有认真更新自己的博客了!正好这段时间在学习异步电机控制,所以把过程中的一些东西写下来,当是回顾也是备忘。本来想是把整个过程的问题和收获都记录下来的(包括硬件设计、mcu控制、算法等),但感觉要好费精力,所以这个异步电机系列就随缘更新吧!自己也是才接触电机控制,有不正确的地方还请多多指教!

        文章中出现的一些名字和概念不理解的话建议可以查查资料哦,博主不可能把所有概念都讲的哦。

二、概述

        我理解的异步电机控制分为开环控制和闭环控制(传统意义上也可以分为标量和矢量),开环一般为v/f控制,较简单,闭环又分为有无传感器。做闭环的话磁链观测器就很重要了,常分为电压型磁链观测器和电流型磁链观测器,而磁链观测器又跟电机参数息息相关,比如定子电阻对电压型磁链观测器的磁链幅值影响很大,所以我们就要对电机参数进行辨识,就引入今天的主题啦!当然你有现成的电机参数就不用辨识了。参数辨识分为在线辨识和离线辨识,今天的主题是离线辨识。

三、原理

        整个工程实现都是参考的一篇硕士论文,最后会把论文贴出来,论文写的很清楚,接下来的原理都是论文上的内容哦。

        论文地址:https://download.csdn.net/download/qq_35141454/87771052

3.1 定子电阻辨识

        原理就是单相直流试验,给电机某一相通一直流,利用电感通直流阻交流的特点,就可以得到电机的等效电路图,求得定子电阻。如下图,

         利用电流闭环先让电流稳定,然后在采集多个电流点,要考虑Mos管导通压降和死区时间的影响,最后通过最小二乘法计算定子电阻,工程实际应用流程如下,

 3.2 转子电阻及漏感辨识

        原理为单相交流实验,其等效电路图与计算原理如下,

 

  3.3 互感辨识

        原理是空载试验,这时同步转速基本等于转子转速,电机转子回路相当于开路了。具体原理如下,

 四、工程代码

        大家应该最关心的是这个了吧,因为理论跟实践是两个东西,有同学肯定感同身受。认真看了代码可能会有疑惑,有些api找不到实现,比如获取三相电流电压,或者某些标志位,又或者一些平台相关的控制代码。这里博主想讲的是这个是一个具体实现思路,目的在于实现,不是要每个细节都体现出来,平台不一样一些东西就不一样(代码是gd32平台),理解了话就很容易了。(为什么会这么说,因为本人在看代码缺少一些具体实现的时候也会很苦恼,哈哈)

        实现的c文件如下,

/*
 * Copyright (c) 2006-2021, RT-Thread Development Team
 *
 * SPDX-License-Identifier: Apache-2.0
 *
 *  @brief: 异步电机参数辨识
 * Change Logs:
 * Date           Author       Notes
 * 2023-04-20     liujian       the first version
 */

#include "param_iden.h"
#include "bus_current.h"
#include <stdlib.h>
#include <string.h>
#include <math.h>

#if 1
// sqrt(3)
#define SQRT3               (1.732051f)
// sqrt(3)/3
#define SQRT3_DIV3           (0.577350f)
// 定子电阻辨识电流环kp
#define RS_KP           (1.0f)
// 定子电阻辨识电流环ki
#define RS_TI           (0.1f)

// Udc值
#define UDC_VALUE               (72)

// 定子电阻辨识
rs_param_t m_rs_param = {0};

// sin0~90度的表
static int16_t g_sin0_90_table[] = {\
    0x0000,0x00C9,0x0192,0x025B,0x0324,0x03ED,0x04B6,0x057F,\
    0x0648,0x0711,0x07D9,0x08A2,0x096A,0x0A33,0x0AFB,0x0BC4,\
    0x0C8C,0x0D54,0x0E1C,0x0EE3,0x0FAB,0x1072,0x113A,0x1201,\
    0x12C8,0x138F,0x1455,0x151C,0x15E2,0x16A8,0x176E,0x1833,\
    0x18F9,0x19BE,0x1A82,0x1B47,0x1C0B,0x1CCF,0x1D93,0x1E57,\
    0x1F1A,0x1FDD,0x209F,0x2161,0x2223,0x22E5,0x23A6,0x2467,\
    0x2528,0x25E8,0x26A8,0x2767,0x2826,0x28E5,0x29A3,0x2A61,\
    0x2B1F,0x2BDC,0x2C99,0x2D55,0x2E11,0x2ECC,0x2F87,0x3041,\
    0x30FB,0x31B5,0x326E,0x3326,0x33DF,0x3496,0x354D,0x3604,\
    0x36BA,0x376F,0x3824,0x38D9,0x398C,0x3A40,0x3AF2,0x3BA5,\
    0x3C56,0x3D07,0x3DB8,0x3E68,0x3F17,0x3FC5,0x4073,0x4121,\
    0x41CE,0x427A,0x4325,0x43D0,0x447A,0x4524,0x45CD,0x4675,\
    0x471C,0x47C3,0x4869,0x490F,0x49B4,0x4A58,0x4AFB,0x4B9D,\
    0x4C3F,0x4CE0,0x4D81,0x4E20,0x4EBF,0x4F5D,0x4FFB,0x5097,\
    0x5133,0x51CE,0x5268,0x5302,0x539B,0x5432,0x54C9,0x5560,\
    0x55F5,0x568A,0x571D,0x57B0,0x5842,0x58D3,0x5964,0x59F3,\
    0x5A82,0x5B0F,0x5B9C,0x5C28,0x5CB3,0x5D3E,0x5DC7,0x5E4F,\
    0x5ED7,0x5F5D,0x5FE3,0x6068,0x60EB,0x616E,0x61F0,0x6271,\
    0x62F1,0x6370,0x63EE,0x646C,0x64E8,0x6563,0x65DD,0x6656,\
    0x66CF,0x6746,0x67BC,0x6832,0x68A6,0x6919,0x698B,0x69FD,\
    0x6A6D,0x6ADC,0x6B4A,0x6BB7,0x6C23,0x6C8E,0x6CF8,0x6D61,\
    0x6DC9,0x6E30,0x6E96,0x6EFB,0x6F5E,0x6FC1,0x7022,0x7083,\
    0x70E2,0x7140,0x719D,0x71F9,0x7254,0x72AE,0x7307,0x735E,\
    0x73B5,0x740A,0x745F,0x74B2,0x7504,0x7555,0x75A5,0x75F3,\
    0x7641,0x768D,0x76D8,0x7722,0x776B,0x77B3,0x77FA,0x783F,\
    0x7884,0x78C7,0x7909,0x794A,0x7989,0x79C8,0x7A05,0x7A41,\
    0x7A7C,0x7AB6,0x7AEE,0x7B26,0x7B5C,0x7B91,0x7BC5,0x7BF8,\
    0x7C29,0x7C59,0x7C88,0x7CB6,0x7CE3,0x7D0E,0x7D39,0x7D62,\
    0x7D89,0x7DB0,0x7DD5,0x7DFA,0x7E1D,0x7E3E,0x7E5F,0x7E7E,\
    0x7E9C,0x7EB9,0x7ED5,0x7EEF,0x7F09,0x7F21,0x7F37,0x7F4D,\
    0x7F61,0x7F74,0x7F86,0x7F97,0x7FA6,0x7FB4,0x7FC1,0x7FCD,\
    0x7FD8,0x7FE1,0x7FE9,0x7FF0,0x7FF5,0x7FF9,0x7FFD,0x7FFE };


/**
 * @brief 等幅值clarke变换
 * @param motor
 * @return
 */
void ea_clarke(ea_clarke_t *model)
{
    model->outa = model->ina;
    model->outb = SQRT3_DIV3 * (model->inb*2 + model->ina);
}

/**
 * @brief pi控制器
 * @param motor
 * @return
 */
void pi_controller(m_pid_t *model)
{
    float tmp = 0;
    float ek = model->target - model->current;

    model->ek_sum = model->ek_sum_last + ek;
    tmp = model->ki * model->ek_sum;
    // 限制积分饱和
    if ((tmp - model->max) > 0.00001f)
        tmp = model->max;
    else if((tmp - model->min) < 0.00001f)
        tmp = model->min;

    model->output = tmp + ek * model->kp;
    // 限制输出
    if ((model->output - model->max) > 0.00001f)
        model->output = model->max;
    else if((model->output - model->min) < 0.00001f)
        model->output = model->min;

    model->ek_sum_last = model->ek_sum;
}

/**
 * @brief svpwm输出通道比较值
 * @param motor
 */
uint16_t ccr1 = 0;
uint16_t ccr2 = 0;
uint16_t ccr3 = 0;
static void motor_svpwm(mrpark_t *target)
{
    // 两个矢量作用时间
    int32_t Tx = 0;
    int32_t Ty = 0;
    // 存放临时比例因子
    float s1 = 0.0f;
    float s2 = 0.0f;
    // ABC三相作用时间
    int32_t pha_time = 0;
    int32_t phb_time = 0;
    int32_t phc_time = 0;
    // 用于计算扇区
    uint8_t A = 0;
    uint8_t B = 0;
    uint8_t C = 0;

    // 计算u1 u2 u3用于判断扇区和统计
    int32_t  U1 = target->u_beta;
    int32_t  U2 = (SQRT3 * target->u_alpha - target->u_beta)/2.0f;
    int32_t  U3 = (-SQRT3 * target->u_alpha - target->u_beta)/2.0f;

    A = U1 > 0? 1:0;
    B = U2 > 0? 1:0;
    C = U3 > 0? 1:0;

    // 计算扇区
    uint8_t N = 4*C + 2*B + A;

    // 将udc等放大32767倍
    float k = SQRT3 * PWM_PERIOD / (UDC_VALUE*32767);

    // 根据扇区各标量计算作用时间
    switch (N)
    {
        case 3: // 扇区1
            Tx = k * U2;// T4
            Ty = k * U1;// T6
        break;

        case 1: // 扇区2
            Tx = -k*U2;// T2
            Ty = -k*U3;// T6
        break;

        case 5: // 扇区3
            Tx = k*U1;// T2
            Ty = k*U3;// T3
        break;

        case 4:
            Tx = -k*U1;// T1
            Ty = -k*U2;// T3
        break;

        case 6: // 扇区5
            Tx = k*U3;// T1
            Ty = k*U2;// T5
        break;

        case 2: // 扇区6
            Tx = -k*U3;// T4
            Ty = -k*U1;// T5
        break;
    }

    // 作用时间大于TS则比例化
    if ((Tx + Ty) > PWM_PERIOD)
    {
        s1 = Tx*1.0/(Tx+Ty);
        s2 = Ty*1.0/(Tx+Ty);
        Tx = s1 * PWM_PERIOD;
        Ty = s2 * PWM_PERIOD;
    }

    // 计算ABC三相作用
    switch( N )
    {
        case 3: // 扇区1
            phc_time = (PWM_PERIOD - Tx - Ty )/2;
            phb_time = phc_time+Ty;
            pha_time = phb_time+Tx;
        break;

        case 1:// 扇区2
            phc_time = (PWM_PERIOD - Tx - Ty )/2;
            pha_time = (phc_time+Ty);
            phb_time = (pha_time+Tx);
        break;

        case 5:// 扇区3
            pha_time = (PWM_PERIOD - Tx - Ty )/2;
            phc_time = pha_time+Ty;
            phb_time = phc_time+Tx;
        break;

        case 4:// 扇区4
            pha_time = (PWM_PERIOD - Tx - Ty )/2;
            phb_time = pha_time+Ty;
            phc_time = phb_time+Tx;
        break;

        case 6:// 扇区5
            phb_time = (PWM_PERIOD - Tx - Ty )/2;
            pha_time = phb_time+Ty;
            phc_time = pha_time+Tx;
        break;

        case 2:// 扇区6
            phb_time = (PWM_PERIOD - Tx - Ty )/2;
            phc_time = phb_time+Ty;
            pha_time = phc_time+Tx;
        break;

        default:
            phb_time = 0;
            phc_time = 0;
            pha_time = 0;
        break;
    }

    // 除以2的原因是PWM为中央对齐模式
    ccr1 = pha_time/2;
    ccr2 = phb_time/2;
    ccr3 = phc_time/2;

    // 设置pwm比较值
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_0, ccr1);
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_1, ccr2);
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_2, ccr3);
}

/**
 * @brief 中值滤波
 * @param motor
 * @return
 */
static float median_filter(float *data, uint16_t len, uint8_t flag)
{
    int i;
    float min, max;
    float sum = 0.0f;
    float median = 0.0f;

    if (NULL == data)
        return 0;

    min = data[0];
    max = data[0];

    for (i=0; i<len; i++)
    {
        sum += data[i];
        if (data[i] - max > 0.000001f) // 找到最大值
        {
            max = data[i];
        }
        else if (data[i] - min < 0.000001f) // 找到最小值
        {
            min = data[i];
        }
    }

    if (!flag)
        median = (sum - min - max)/(len-2); // 去掉最大和最小求平均值
    else
        median = sum / len; // 直接求平均值
  //  rt_kprintf("min=%f, max=%f, median=%f\n", min, max, median);
    return median;
}

/**
 * @brief 最小二乘法
 *        y =ax+b  a=x
 * @param motor
 * @return
 */
static float least_square(float *data, uint16_t len)
{

}

/**
 * @brief 定子参数辨识
 * @param motor
 * @return
 */
/* 目标电流个数 */
#define TARGET_CURRENT_CNT      (15)
/* clarke变换 */
ea_clarke_t clarke;
/* 缓存vbus */
float vbus_tmp = 0.0f;
/* 目标电流 */
static float g_current_target_table[TARGET_CURRENT_CNT] = {1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8
                                         ,3.0,3.2,3.4,3.6,3.8};
void current_loop(void)
{
    mrpark_t park;

    // 获取三相电流
    get_abc_current((float *)&clarke);
    // 获取母线电压
    vbus_tmp = get_vm_bus_voltage();
    // 等幅值clarke变换
    ea_clarke(&clarke);
    // pi调节
    m_rs_param.pid.current = -clarke.outa;
    pi_controller(&m_rs_param.pid);

    park.u_alpha = m_rs_param.pid.output*32767;
    park.u_beta = 0;
    // svpwm输出
    motor_svpwm(&park);
}

/**
 * @brief 定子参数辨识
 *          为adc采集周期调用
 *          step1:等待目标电流稳定
 *          step2:计算单个目标电流均值
 *          step3:最小二乘法计算定子电阻
 * @param None
 * @return 0
 */
int iden_stator_resistance(void)
{
    uint32_t i;
    float vbus = 0.0f;

    float tmp[TARGET_CURRENT_CNT] = {0};

    // 电流闭环调节
    current_loop();

    switch (m_rs_param.state) {
    case RS_LEVEL_STANDBY:  // 等待目标电流稳定下来

        m_rs_param.cnt++;

        if (m_rs_param.cnt > PWM_HZ/10) // 100ms后开始计算
        {
            m_rs_param.cnt = 0;
        //    m_rs_param.mcnt = 0;
            m_rs_param.state = RS_LEVEL_CALC1;
        }
        break;

    case RS_LEVEL_CALC1:   // 采样100个电压电流点,进行数字滤波

        m_rs_param.isam1[m_rs_param.cnt] = -clarke.ina;
        // 减去死区时间 deltaT/(2Ts)*udc的压降 减去Mos管导通压降
        vbus = vbus_tmp * ((ccr1*2 - ccr2 - ccr3)/3.0/(PWM_PERIOD/2))
                                - (PWM_HZ*0.17/2000000*vbus_tmp); // 死区时间0.17us
                                - 0.0093*(-clarke.ina); // 导通电阻0.0093

        m_rs_param.usam1[m_rs_param.cnt] = vbus;
        m_rs_param.cnt++;

        if (m_rs_param.cnt >= RS_SAMP_CNT)
        {
            // 中值滤波
            m_rs_param.isam2[m_rs_param.mcnt] = median_filter(m_rs_param.isam1, RS_SAMP_CNT, 0);
            m_rs_param.usam2[m_rs_param.mcnt] = median_filter(m_rs_param.usam1, RS_SAMP_CNT, 0);
            m_rs_param.mcnt++;

            // 目标电流集合采集计算完毕
            if (m_rs_param.mcnt >= TARGET_CURRENT_CNT)
            {
                m_rs_param.mcnt = 0;
                m_rs_param.state = RS_LEVEL_CALC2;
            }
            else // 计算下一个目标电流
            {
                m_rs_param.pid.target = m_rs_param.itable[m_rs_param.mcnt];
                m_rs_param.state = RS_LEVEL_STANDBY;
                m_rs_param.cnt = 0;
            }
        }
        break;
    case RS_LEVEL_CALC2:

        // 中值滤波
        m_rs_param.xmedian = median_filter(m_rs_param.isam2, TARGET_CURRENT_CNT, 1); // 计算x平均
        m_rs_param.ymedian = median_filter(m_rs_param.usam2, TARGET_CURRENT_CNT, 1);// 计算y平均

        for (i=0; i<TARGET_CURRENT_CNT; i++)
            tmp[i] = m_rs_param.isam2[i] * m_rs_param.isam2[i];
        m_rs_param.x2median = median_filter(tmp, TARGET_CURRENT_CNT, 1); // 计算x^2平均

        for (i=0; i<TARGET_CURRENT_CNT; i++)
            tmp[i] = m_rs_param.usam2[i] * m_rs_param.isam2[i];
        m_rs_param.xymedian = median_filter(tmp, TARGET_CURRENT_CNT, 1);// 计算xy平均


        // 最小二乘法求得定子电阻,电压与电流关系如下 U=ai+b,a为定子电阻
//        m_rs_param.rs_res = ( (TARGET_CURRENT_CNT)*m_rs_param.xymedian- m_rs_param.xmedian* m_rs_param.ymedian)
//                            / ((TARGET_CURRENT_CNT)*m_rs_param.x2median - m_rs_param.xmedian * m_rs_param.xmedian);
        m_rs_param.rs_res = m_rs_param.xymedian / m_rs_param.x2median;
        m_rs_param.bres = m_rs_param.ymedian - (m_rs_param.xmedian*m_rs_param.rs_res);

        // 关闭三项桥输出
        timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_0, 0);
        timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_1, 0);
        timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_2, 0);

        // 停止计算
        m_rs_param.state = RS_LEVEL_NONE;
        m_rs_param.rs_start_flag = 0;
        break;
    default:break;
    }

    return 0;
}



/**
 * @brief 定子参数辨识
 * @param motor
 * @return
 */
//int rs_start(int argc, char **argv)
int rs_start(void)
{
    uint32_t cnt = 0;
    float current[3] = {0};
    float  target = 0;
/*
    if (argc < 2)
    {
        rt_kprintf("parameter < 2!\n");
        return -1;
    }

    target = atof(argv[1]);
*/
    memset(&m_rs_param, 0, sizeof(rs_param_t));

    m_rs_param.pid.kp = RS_KP;
    m_rs_param.pid.ki = RS_TI;
    m_rs_param.pid.min = 0.01;
    m_rs_param.pid.max = 3.0; // 定子电阻很小,目标不能给很大,否则有烧坏电机和板子的风险

    m_rs_param.i1_target = 1.2f;
    m_rs_param.i2_target = target;//3.0f;

    // 关闭三项桥
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_0, 0);
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_1, 0);
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_2, 0);
    // 等待电流稳定后计算偏置电流
    rt_thread_mdelay(20);
    // 计算偏置电流
    calc_phase_offset();

    // 各个参数初始化
    m_rs_param.state = RS_LEVEL_STANDBY;
    m_rs_param.cnt = 0;
    m_rs_param.mcnt = 0;
    m_rs_param.itable = g_current_target_table; // 这个表的参考电流能给太大,否则烧坏电机和板子
    m_rs_param.pid.target = g_current_target_table[0];//m_rs_param.i2_target;
    m_rs_param.rs_start_flag = 1;
#if 1
    rt_thread_mdelay(2000);
    rt_kprintf("Rs=%f\n", m_rs_param.rs_res);
#else

     // 等待参数计算完毕
    rt_thread_mdelay(2000);

    while(1)
    {
     //   get_abc_current(current);
      //  rt_kprintf("%d %f %f %f\n", cnt, current[0], current[1], current[2]);
      //  rt_kprintf("%d %f %f %f %f %d %d %d\n", cnt, clarke.ina, clarke.inb, clarke.inc,
      //          m_rs_param.pid.output, ccr1, ccr2,ccr3);
        rt_kprintf("%d %f %f %f\n", cnt, m_rs_param.isam2[cnt], m_rs_param.usam2[cnt], m_rs_param.usam2[cnt]/m_rs_param.isam2[cnt]);
        cnt++;
        if (cnt >= TARGET_CURRENT_CNT)
        {
            rt_kprintf("clac: %f %f %f %f %f %f\n", m_rs_param.x2median, m_rs_param.xymedian,
                    m_rs_param.xmedian, m_rs_param.ymedian, m_rs_param.bres,m_rs_param.rs_res);
            break;
        }
        rt_thread_mdelay(10);
    }
#endif
    return 0;
}
//MSH_CMD_EXPORT(rs_start, rs start);






rr_param_t g_rr_param = {0};

/**
 * @brief 获取当前角度下的sin cos
 * @param motor
 */

void get_sin_cos(theta_t *dtc)
{
    uint8_t    sector   = 0;
    uint16_t   index  = 0;

    index = (dtc->theta) >> 6;            //< 0~65535等比例缩小64倍,范围0~1024

    //根据扇区计算sin和cos值
    sector = (index & 0x300) >> 8;            //0:0~255 1:256~511 2:512~767 3:767~1023
    switch(sector)
    {
        case 0:// 0~90度
            dtc->sin = (g_sin0_90_table[(uint8_t)(index)]);
            dtc->cos = (g_sin0_90_table[(uint8_t)(0xFFu - (uint8_t)(index))]);
        break;

        case 1:// 90~180度
            dtc->sin = (g_sin0_90_table[(uint8_t)(0xFFu - (uint8_t)(index))]);
            dtc->cos = (-g_sin0_90_table[(uint8_t)(index)]);
        break;

        case 2://  180~270度
            dtc->sin = ( -g_sin0_90_table[(uint8_t)(index)]);
            dtc->cos = (-g_sin0_90_table[(uint8_t)(0xFFu - (uint8_t)(index))]);
        break;

        case 3:// 270~360度
            dtc->sin = (-g_sin0_90_table[(uint8_t)(0xFFu - (uint8_t)(index))]);
            dtc->cos = (g_sin0_90_table[(uint8_t)(index)]);
        break;
    }
}

static complex result = {0};
/**
 * @brief 离散傅里叶变换
 * @param
 */
void dft_calc(float *input, int k, uint16_t num)
{
    int n = 0;
    complex sum;
    complex point;//[RR_SAM_CNT];

    sum.real = 0;
    sum.imag = 0;
    uint16_t N = num;

    for(n=0; n<N; n++)
    {
//        point[n].real = cos(2*M_PI/N*k*n)*input[n];  //运用欧拉公式把复数拆分成实部和虚部
//        point[n].imag = -sin(2*M_PI/N*k*n)*input[n];
//        sum.real += point[n].real;//对每一个点的实部求和
//        sum.imag += point[n].imag;
        point.real = cos(2*M_PI/N*k*n)*input[n];  //运用欧拉公式把复数拆分成实部和虚部
        point.imag = -sin(2*M_PI/N*k*n)*input[n];

        sum.real += point.real;//对每一个点的实部求和
        sum.imag += point.imag;
        //对每一个点的虚部求和
    }

    result.real = sum.real;
    result.imag = sum.imag;
}


/**
 * @brief 转子电阻和漏感参数辨识
 *          adc中断调用
 * @param None
 * @return 0
 */
int irr_current_intr_handler(void)
{
    float i_tmp[3] = {0};
    float vbus_tmp = 0.0f;
    static float last_i = 0.0f;
    float cur = 0.0f;

    // 获取三相电流
    get_abc_current(i_tmp);
    // 获取母线电压
    vbus_tmp = get_vm_bus_voltage();

    // 低通滤波  这样会导致相位后移
  //  cur = ((i_tmp[0] / 10000.0 + 0.0003 * last_i) / (0.0004));
   // last_i = cur;

    g_rr_param.isam[g_rr_param.cnt] = -i_tmp[0];//cur;
    g_rr_param.usam[g_rr_param.cnt] = vbus_tmp * (ccr1 - ccr2)/(PWM_PERIOD/2.0) //Va-Vb - 死区电压-Mos管压降
                                        - (PWM_HZ*0.17/2000000*vbus_tmp)
                                        - 0.0093*(i_tmp[0]);
    g_rr_param.cnt++;

    // 采集一周期数据完成
    if (g_rr_param.cnt >= RR_SAM_CNT)
    {
        g_rr_param.rr_start_flag = 0;
        g_rr_param.rr_start_flag1 = 0;
        g_rr_param.cnt = 0;
        timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_0, 0);
        timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_1, 0);
        timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_2, 0);

        // 电流傅里叶变换  (2/N)*(求和x(n)cos(2*pi*n*k/N)-j求和x(n)sin(2*pi*n*k/N))
        dft_calc(g_rr_param.isam, 1, RR_SAM_CNT); //k=1,求的是基波频率的傅里叶变换
        g_rr_param.i_dft.real = result.real*2/RR_SAM_CNT;
        g_rr_param.i_dft.imag = result.imag*2/RR_SAM_CNT;
        // 求电流基波分量幅值
        g_rr_param.i_amp = sqrt(g_rr_param.i_dft.real*g_rr_param.i_dft.real
                                + g_rr_param.i_dft.imag*g_rr_param.i_dft.imag);

        // 电压傅里叶变换
        dft_calc(g_rr_param.usam, 1, RR_SAM_CNT); //k=1,求的是基波频率的傅里叶变换
        g_rr_param.u_dft.real = result.real*2/RR_SAM_CNT;
        g_rr_param.u_dft.imag = result.imag*2/RR_SAM_CNT;
        // 求电压基波分量幅值
        g_rr_param.u_amp = sqrt(g_rr_param.u_dft.real*g_rr_param.u_dft.real
                                + g_rr_param.u_dft.imag*g_rr_param.u_dft.imag);

        // 求功率因数角cos值( Ua*Ia + Ub*Ib)/u_amp*i_amp
        g_rr_param.acos = ((g_rr_param.i_dft.real*g_rr_param.u_dft.real)
                         +(g_rr_param.i_dft.imag*g_rr_param.u_dft.imag))
                                 /(g_rr_param.u_amp * g_rr_param.i_amp);
        // 求功率因数角sin值( Ub*Ia - Ua*Ib)/u_amp*i_amp
        g_rr_param.asin = ((g_rr_param.u_dft.imag*g_rr_param.i_dft.real)
                         - (g_rr_param.u_dft.real*g_rr_param.i_dft.imag))
                                 /(g_rr_param.u_amp * g_rr_param.i_amp);

        // 求得转子电阻 Rr=(Vab/Ia)cos*(2/3)-Rs
        g_rr_param.rr = (g_rr_param.u_amp/g_rr_param.i_amp)*g_rr_param.acos*2/3.0;
        // 求得转子漏感和定子漏感 lls=llr=(Vab/Ia)sin/3/2*pi*f
        g_rr_param.lls = (g_rr_param.u_amp/g_rr_param.i_amp)*g_rr_param.asin/3.0/(2*M_PI*BASE_VOL_FREQ);
        g_rr_param.llr = g_rr_param.lls;
    }
}

/**
 * @brief 转子电阻和漏感参数辨识
 *          pwm更新中断调用
 * @param None
 * @return 0
 */
int irr_pwm_intr_handler(void)
{
    mrpark_t park;
    // 每个pwm周期转了多少转
    float cycle_freq = ((float)(BASE_VOL_FREQ)/(float)PWM_HZ);

    g_rr_param.theta += round(cycle_freq * 65535); // 放大65535倍,每次加freq_ * 65535

    get_sin_cos((theta_t *)&g_rr_param);

    park.u_alpha = BASE_VAL_AMP * g_rr_param.cos;
    park.u_beta = 0;
    // svpwm输出
    motor_svpwm(&park);
}

//int rr_start(int argc, char **argv)
int rr_start(void)
{
    int cnt = 0;

    memset(&g_rr_param, 0, sizeof(rr_param_t));

    // 关闭三相桥
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_0, 0);
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_1, 0);
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_2, 0);
    // 等待电流稳定后计算偏置电流
    rt_thread_mdelay(20);
    // 计算偏置电流
    calc_phase_offset();

    // 开启pwm
    g_rr_param.rr_start_flag = 1;
    // 等待输出稳定
    rt_thread_mdelay(200);
    // 开启adc采集
    g_rr_param.rr_start_flag1 = 1;
    // 等待数据结果
    rt_thread_mdelay(2000);
#if 1
  //  g_rr_param.rr = g_rr_param.rr - m_rs_param.rs_res;
    rt_kprintf("Rr=%f llr=%f lls=%f\n", g_rr_param.rr, g_rr_param.llr,  g_rr_param.llr);
#else
    while (1)
    {
        rt_kprintf("%d %f %f %d %d\n", cnt, g_rr_param.isam[cnt], g_rr_param.usam[cnt], ccr1, ccr2);
        cnt++;
        if (cnt >= RR_SAM_CNT)
        {

            rt_kprintf("total:%f %f %f %f %f %f %f\n", g_rr_param.i_dft.real, g_rr_param.i_dft.imag
                    , g_rr_param.i_amp, g_rr_param.u_dft.real,g_rr_param.u_dft.imag,g_rr_param.u_amp
                    ,g_rr_param.acos);
            rt_kprintf("Rr=%f llr=%f\n", g_rr_param.rr, g_rr_param.llr);
            break;
        }
        rt_thread_mdelay(10);
    }
#endif
}
//MSH_CMD_EXPORT(rr_start, rr start);



lm_param_t g_lm_param = {0};

/**
 * @brief 互感参数辨识
 *          pwm更新中断调用
 * @param None
 * @return 0
 */
int ilm_pwm_intr_handler(void)
{
    mrpark_t park;
    // 每个pwm周期转了多少转
    float cycle_freq = ((float)(ILM_BASE_VOL_FREQ)/(float)PWM_HZ);

    g_lm_param.theta += round(cycle_freq * 65535); // 放大65535倍,每次加freq_ * 65535

    get_sin_cos((theta_t *)&g_lm_param);

    park.u_alpha = ILM_BASE_VAL_AMP * g_lm_param.cos;
    park.u_beta = ILM_BASE_VAL_AMP * g_lm_param.sin;
    // svpwm输出
    motor_svpwm(&park);
}

/**
 * @brief 互感参数辨识
 *          adc中断调用
 * @param None
 * @return 0
 */
int ilm_current_intr_handler(void)
{
    float i_tmp[3] = {0};
    float vbus_tmp = 0.0f;
    static float last_i = 0.0f;
    float cur = 0.0f;

    // 获取三相电流
    get_abc_current(i_tmp);
    // 获取母线电压
    vbus_tmp = get_vm_bus_voltage();

    // 低通滤波  这样会导致相位后移
  //  cur = ((i_tmp[0] / 10000.0 + 0.0003 * last_i) / (0.0004));
   // last_i = cur;

    g_lm_param.isam[g_lm_param.cnt] = -i_tmp[0];//cur;
    g_lm_param.usam[g_lm_param.cnt] = vbus_tmp * ((ccr1*2 - ccr2 - ccr3)/3.0/(PWM_PERIOD/2)); //Va - 死区电压-Mos管压降
                                      //  - (PWM_FREQUENCY*0.17/2000000*vbus_tmp)
                                       // - 0.0093*(i_tmp[0]);;
    g_lm_param.cnt++;

    // 采集一周期数据完成
    if (g_lm_param.cnt >= LM_SAM_CNT)
    {
        g_lm_param.lm_start_flag = 0;
        g_lm_param.lm_start_flag1 = 0;
        g_lm_param.cnt = 0;
        timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_0, 0);
        timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_1, 0);
        timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_2, 0);

        // 电流傅里叶变换  (2/N)*(求和x(n)cos(2*pi*n*k/N)-j求和x(n)sin(2*pi*n*k/N))
        dft_calc(g_lm_param.isam, 1, LM_SAM_CNT); //k=1,求的是基波频率的傅里叶变换
        g_lm_param.i_dft.real = result.real*2.0/LM_SAM_CNT;
        g_lm_param.i_dft.imag = result.imag*2.0/LM_SAM_CNT;
        // 求电流基波分量幅值
        g_lm_param.i_amp = sqrt(g_lm_param.i_dft.real*g_lm_param.i_dft.real
                                + g_lm_param.i_dft.imag*g_lm_param.i_dft.imag);

        // 电压傅里叶变换
        dft_calc(g_lm_param.usam, 1, LM_SAM_CNT); //k=1,求的是基波频率的傅里叶变换
        g_lm_param.u_dft.real = result.real*2.0/LM_SAM_CNT;
        g_lm_param.u_dft.imag = result.imag*2.0/LM_SAM_CNT;
        // 求电压基波分量幅值
        g_lm_param.u_amp = sqrt(g_lm_param.u_dft.real*g_lm_param.u_dft.real
                                + g_lm_param.u_dft.imag*g_lm_param.u_dft.imag);

        // 求功率因数角cos值( Ua*Ia + Ub*Ib)/u_amp*i_amp
//        g_lm_param.acos = ((g_lm_param.i_dft.real*g_lm_param.u_dft.real)
//                         +(g_lm_param.i_dft.imag*g_lm_param.u_dft.imag))
//                                 /(g_lm_param.u_amp * g_lm_param.i_amp);
        // 求功率因数角sin值( Ub*Ia - Ua*Ib)/u_amp*i_amp
        g_lm_param.asin = ((g_lm_param.u_dft.imag*g_lm_param.i_dft.real)
                         - (g_lm_param.u_dft.real*g_lm_param.i_dft.imag))
                                 /(g_lm_param.u_amp * g_lm_param.i_amp);

        // 求得转子漏感和定子漏感 lls=llr=(Vab/Ia)sin/3/2*pi*f
        g_lm_param.lm = (g_lm_param.u_amp/g_lm_param.i_amp)*g_lm_param.asin/(2*M_PI*ILM_BASE_VOL_FREQ);
    }
}

//int lm_start(int argc, char **argv)
int lm_start(void)
{
    int cnt = 0;

    memset(&g_lm_param, 0, sizeof(lm_param_t));

    // 关闭三相桥
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_0, 0);
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_1, 0);
    timer_channel_output_pulse_value_config(TIMER7, TIMER_CH_2, 0);
    // 等待电流稳定后计算偏置电流
    rt_thread_mdelay(20);
    // 计算偏置电流
    calc_phase_offset();

    // 开启pwm
    g_lm_param.lm_start_flag = 1;
    // 等待输出稳定
    rt_thread_mdelay(3000);
    // 开启adc采集
    g_lm_param.lm_start_flag1 = 1;
    // 等待数据结果
    rt_thread_mdelay(2000);
#if 1
    rt_kprintf("Lm=%f\n", g_lm_param.lm);
#else
    while (1)
    {
        rt_kprintf("%d %f %f %d %d\n", cnt, g_lm_param.isam[cnt], g_lm_param.usam[cnt], ccr1, ccr2);
        cnt++;
        if (cnt >= LM_SAM_CNT)
        {

            rt_kprintf("total:%f %f %f %f %f %f %f\n", g_lm_param.i_dft.real, g_lm_param.i_dft.imag
                    , g_lm_param.i_amp, g_lm_param.u_dft.real,g_lm_param.u_dft.imag,g_lm_param.u_amp
                    ,g_lm_param.asin);
            rt_kprintf("Lm=%f\n", g_lm_param.lm);
            break;
        }
        rt_thread_mdelay(10);
    }
#endif
}
//MSH_CMD_EXPORT(lm_start, lm start);
#if 1



motor_iden_param_t param = {0};

motor_iden_param_t *get_motor_iden_param(void)
{
    return &param;
}

int start_iden(int argc, char **argv)
{
    int i;
    int cnt = 0;
    float rs_sum = 0.0f;
    float rr_sum = 0.0f;
    float llr_sum = 0.0f;
    float lls_sum = 0.0f;
    float lm_sum = 0.0f;

    if (argc < 2)
    {
        rt_kprintf("argc <2!\n");
        return -1;
    }

    cnt = atoi(argv[1]);

    for (i=0; i<cnt; i++)
    {
        rs_start();
        rt_thread_mdelay(100);

        rr_start();
        g_rr_param.rr = g_rr_param.rr - m_rs_param.rs_res;

        rt_thread_mdelay(100);
        lm_start();
        g_lm_param.lm = g_lm_param.lm - g_rr_param.llr;


      //  rt_kprintf("*********************************iden total*********************************\n");
        rt_kprintf("****(%d):     Rs=%f,Rr=%f,Llr=Lls=%f,Lm=%f\n",i+1, m_rs_param.rs_res, g_rr_param.rr
                ,g_rr_param.llr, g_lm_param.lm);
       // rt_kprintf("************************************end*************************************\n");

        rs_sum += m_rs_param.rs_res;
        rr_sum += g_rr_param.rr;
        llr_sum += g_rr_param.llr;
        lls_sum = llr_sum;
        lm_sum += g_lm_param.lm;

        rt_thread_mdelay(100);
    }

   // rt_base_t level = rt_hw_interrupt_disable();
    param.calc_rs = rs_sum/cnt;
    param.calc_rr = rr_sum/cnt;
    param.calc_lm = lm_sum/cnt;
    param.calc_ls = lls_sum/cnt + param.calc_lm;
    param.calc_lr = llr_sum/cnt + param.calc_lm;
    param.calc_sigma = 1-(param.calc_lm*param.calc_lm)/((param.calc_ls)*(param.calc_lr));
  //  rt_hw_interrupt_enable(level);

    rt_kprintf("**************iden param***********\n");
    rt_kprintf("*       Rs:  %f             *\n", param.calc_rs);
    rt_kprintf("*       Rr:  %f             *\n", param.calc_rr);
    rt_kprintf("*       ls:  %f            *\n", param.calc_ls);
    rt_kprintf("*       lr:  %f            *\n", param.calc_lr);
    rt_kprintf("*       Lm:  %f             *\n", param.calc_lm);
    rt_kprintf("*       sig:  %f            *\n", param.calc_sigma);
    rt_kprintf("*       llr:  %f            *\n", llr_sum/cnt);
    rt_kprintf("**************iden end************\n");
}
MSH_CMD_EXPORT(start_iden, iden param);
#endif
#endif

头文件如下,

/*
 * Copyright (c) 2006-2021, RT-Thread Development Team
 *
 * SPDX-License-Identifier: Apache-2.0
 *
 * Change Logs:
 * Date           Author       Notes
 * 2023-04-20     lj       the first version
 */
#ifndef APPLICATIONS_THREE_PHASE_MOTOR_PARAM_IDEN_H_
#define APPLICATIONS_THREE_PHASE_MOTOR_PARAM_IDEN_H_

#include <stdint.h>

// 设置pwm频率
#define PWM_HZ           (10000)
// 频率
#define CALC_FREQ       PWM_HZ
// pwm周期 时钟/频率
#define PWM_PERIOD              (uint16_t)(200000000/PWM_HZ)
// 基波电压频率100hz
#define BASE_VOL_FREQ   (200)
// 基波电压幅值
#define BASE_VAL_AMP    (2.0f)

// 基波电压频率200hz
#define ILM_BASE_VOL_FREQ   (200)
// 基波电压幅值
#define ILM_BASE_VAL_AMP    (7.5f)

// RS采样次数
#define RS_SAMP_CNT     (100)
// 单相交流采样点数 PWM_HZ/BASE_VOL_FREQ
#define RR_SAM_CNT      (PWM_HZ/BASE_VOL_FREQ) //(100)
// 空载试验采样点数 PWM_HZ/ILM_BASE_VOL_FREQ
#define LM_SAM_CNT      (PWM_HZ/ILM_BASE_VOL_FREQ)//(100)

typedef enum {
    RS_LEVEL_NONE,
    RS_LEVEL_STANDBY, // 电流稳定
    RS_LEVEL_CALC1,   // 计算第一次电流
    RS_LEVEL_CALC2,   // 计算第二次电流
}rs_state_t;

/**
 * @brief pwm通道比较值
 */
typedef struct{
    uint16_t ccr1;  //pwm通道1比较值
    uint16_t ccr2;  //pwm通道2比较值
    uint16_t ccr3;  //pwm通道3比较值
}ccr_t;

typedef struct {
    uint16_t theta;
    float sin;
    float cos;
}theta_t;

/**
 * @brief 反park变化计算得到的alpha_beta
 */
typedef struct {
    int32_t u_alpha;
    int32_t  u_beta;
}mrpark_t;


typedef struct {
    float ina;
    float inb;
    float inc;
    float outa;
    float outb;
}ea_clarke_t;

/**
 * @brief pid控制
 */
typedef struct {
    float current;    //  输入,当前
    float target;    //  输入,目标
    float output;     //  输出,d轴分量目标电流
    float kp;       //  比例
    float ki;       //  积分
    float min;      //  最小值
    float max;      //  最大值
    float ek_sum;   //  误差积分
    float ek_sum_last;   //  误差积分
}m_pid_t;


typedef struct {
    float i1_target;
    float i2_target;
    m_pid_t pid;
    float *itable;
    float isam1[RS_SAMP_CNT];
    float usam1[RS_SAMP_CNT];
    float isam2[RS_SAMP_CNT];
    float usam2[RS_SAMP_CNT];
    float xmedian;
    float ymedian;
    float x2median;
    float xymedian;
   // float imedian[10];
   // float umedian[10];
    uint8_t rs_start_flag;
    uint8_t state;
    uint16_t cnt;    // 记数节拍
    uint16_t mcnt;  // calc2记数节拍
    float rs_res;   //定子电阻
    float bres;   //定子电阻
}rs_param_t;

//复数结构体,用于实现傅里叶运算
typedef struct
{
    double real,imag;
}complex;

typedef struct {
    uint16_t theta;
    float sin;
    float cos;
    float isam[RR_SAM_CNT];
    float usam[RR_SAM_CNT];
    complex i_dft;
    complex u_dft;
    float i_amp;    // 电流幅值
    float u_amp;    // 电流幅值
    float acos;     // 功率因数角cos值
    float asin;     // 功率因数角sin值
    float lls;
    float llr;
    float rr;
    uint16_t cnt; //记数节拍
    uint8_t rr_start_flag;
    uint8_t rr_start_flag1;
}rr_param_t;

typedef struct {
    uint16_t theta;
    float sin;
    float cos;
    float isam[LM_SAM_CNT];
    float usam[LM_SAM_CNT];
    complex i_dft;
    complex u_dft;
    float i_amp;    // 电流幅值
    float u_amp;    // 电流幅值
    float acos;     // 功率因数角cos值
    float asin;     // 功率因数角sin值
    float lm;
    uint16_t cnt; //记数节拍
    uint8_t lm_start_flag;
    uint8_t lm_start_flag1;
}lm_param_t;

typedef struct {
    float calc_rs;
    float calc_rr;
    float calc_ls;
    float calc_lr;
    float calc_lm;
    float calc_sigma;
}motor_iden_param_t;

extern rs_param_t m_rs_param;
extern rr_param_t g_rr_param;
extern lm_param_t g_lm_param;
motor_iden_param_t *get_motor_iden_param(void);

#endif /* APPLICATIONS_THREE_PHASE_MOTOR_PARAM_IDEN_H_ */

  • 6
    点赞
  • 52
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值