如何辨识RM云台黑箱系统？怎么建模？

基础知识

系统辨识及黑箱方法

系统辨识是通过考察系统的输入、输出及其动态过程，来定量或定性地认识系统的功能特性、行为方式，以及探索其内部结构和机理的一种控制论认识方法。
通过辨识建立数学模型的目的是估计表征系统行为的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入和输出预测系统输出的未来演变，以及设计控制器。
黑箱方法，指内部构造还不清楚，由于条件的限制，只能通过外部观测和试验去认识其功能和特性的系统。人们把外部对黑箱的影响称为黑箱的输入，把黑箱对外部的反应称为黑箱的输出。
状态空间方程

系统的动态特性是用由状态变量构成的一阶微分方程组来描述的。它能反映系统的全部独立变量的变化，从而能同时确定系统的全部内部运动状态，而且还可以方便地处理初始条件。
状态空间表达式： 由状态方程和输出方程构成，在状态空间中对控制系统作完整表述的公式。
系统的状态方程为：





其中，A 为 n×n 的常系数矩阵，称作系统矩阵 ； B 为 n×r 的常系数矩阵，称作控制矩阵。 A 与 B 都由系统本身的参数决定。 u 是输入信号， x 是系统状态向量。
系统的输出方程：




其中，C 为 m×n 的常系数矩阵，称为输出矩阵，它表达了输出变量与状态变量之间的关系； D 为 m×r 的常系数矩阵，称为直接转移矩阵，它表示输入变量通过矩阵 D 直接转移到输出。在大多数实际系统中， D=0 。 y 是输出信号。
黑箱系统辨识

基本流程如下图所示






激励信号
利用MATLAB产生频率从1Hz到500Hz变化的扫频信号。由于频率范围较大，且高频部分云台振幅变化较小，为了减少数据量，采用类指数形式的变化趋势。为了保证信号不会发生跳变，在每次信号频率发生变化时都需要从零相位开始，因此周期需要是整数。


clear; clc; close all;


F = ([1:0.5:22. 24:2:40, 50:10:120, 200, 250, 333, 500]);
T = round(1000./F);%周期
F = 1000./T;%频率


plot(F);
ylabel('F/Hz');


这里产生的频率和对应周期是1Hz到500Hz，但经过实测发现用不到过于高频的信号，实际使用的是1Hz到40Hz的信号并达到了更好的效果，代码做一些微调整就行。


云台YAW轴
根据经验和时域的响应情况，调节PID参数构建速度环和位置环使云台可以稳定的跟随目标角速度。


将上一步Matlab生成的F和T导入程序，并产生激励信号驱动电机。这里由于激励信号数据量很大，正弦信号在定时器里产生并驱动。这里我使用了两种方案：使用硬定时器或使用基于FreeRTOS的软定时器。


硬定时器方案：在CubeMX配置好定时器频率为500Hz，将信号生成函数扔进回调函数。这里使用了TIM8。
在tim.c文件中，定时器设置如下：


void MX_TIM8_Init(void)
{
  TIM_ClockConfigTypeDef sClockSourceConfig = {0};
  TIM_MasterConfigTypeDef sMasterConfig = {0};


  htim8.Instance = TIM8;
  htim8.Init.Prescaler = 16800-1;
  htim8.Init.CounterMode = TIM_COUNTERMODE_UP;
  htim8.Init.Period = 20-1;
  htim8.Init.ClockDivision = TIM_CLOCKDIVISION_DIV1;
  htim8.Init.RepetitionCounter = 0;
  htim8.Init.AutoReloadPreload = TIM_AUTORELOAD_PRELOAD_DISABLE;
  if (HAL_TIM_Base_Init(&htim8) != HAL_OK)
  {
    Error_Handler();
  }
  sClockSourceConfig.ClockSource = TIM_CLOCKSOURCE_INTERNAL;
  if (HAL_TIM_ConfigClockSource(&htim8, &sClockSourceConfig) != HAL_OK)
  {
    Error_Handler();
  }
  sMasterConfig.MasterOutputTrigger = TIM_TRGO_RESET;
  sMasterConfig.MasterSlaveMode = TIM_MASTERSLAVEMODE_DISABLE;
  if (HAL_TIMEx_MasterConfigSynchronization(&htim8, &sMasterConfig) != HAL_OK)
  {
    Error_Handler();
  }


}


TIM8设置128000000 / (12800 * 20) = 500 Hz，也就是定时器的频率是500Hz。
回调函数在stm32f4xx_it.c文件中：


#define pi acos(-1)


//float t[64] = {1000, 667, 500, 400, 333, 286, 250, 222, 200, 182, 167, 154, 143, 133, 125, 118, 111, 105, 100, 95, 91, 87, 83, 80, 77, 74, 71, 69, 67, 65, 63, 61, 59, 57, 56, 54, 53, 51, 50, 49, 48, 47, 45, 42, 38, 36, 33, 31, 29, 28, 26, 25, 20, 17, 14, 13, 11, 10, 9, 8, 5, 4, 3, 2};
//float f[64] = {1.000000, 1.499250, 2.000000, 2.500000, 3.003003, 3.496503, 4.000000, 4.504505, 5.000000, 5.494505, 5.988024, 6.493506, 6.993007, 7.518797, 8.000000, 8.474576, 9.009009, 9.523810, 10.000000, 10.526316, 10.989011, 11.494253, 12.048193, 12.500000, 12.987013, 13.513514, 14.084507, 14.492754, 14.925373, 15.384615, 15.873016, 16.393443, 16.949153, 17.543860, 17.857143, 18.518519, 18.867925, 19.607843, 20.000000, 20.408163, 20.833333, 21.276596, 22.222222, 23.809524, 26.315789, 27.777778, 30.303030, 32.258065, 34.482759, 35.714286, 38.461538, 40.000000, 50.000000, 58.823529, 71.428571, 76.923077, 90.909091, 100.000000, 111.111111, 125.000000, 200.000000, 250.000000, 333.333333, 500.000000};


float tt[64];
int sin_signal_switch;
int i = 0;
float time = 0;
float sin_signal;


void TIM8_UP_TIM13_IRQHandler(void)
{
  /* USER CODE BEGIN TIM8_UP_TIM13_IRQn 0 */
        while(time * 1000 > tt && i < 63)
        {
                i++;
        }
        if(i == 64)
        {
                sin_signal = 0;
                time = 0;
                sin_signal_switch = 0;
        }
        sin_signal = sin(2 * pi * f * time);//正弦信号
        //sin_signal = sin(2 * pi * f * time) > 0?0:1;//阶跃信号
        time += 0.002;
  /* USER CODE END TIM8_UP_TIM13_IRQn 0 */
  HAL_TIM_IRQHandler(&htim8);
  /* USER CODE BEGIN TIM8_UP_TIM13_IRQn 1 */


  /* USER CODE END TIM8_UP_TIM13_IRQn 1 */
}


软定时器方案：移植FreeRTOS后，创建一个频率为500Hz的软定时器并执行信号生成函数。


xTaskCreate((TaskFunction_t )signal_generate,            
            (const char*    )"signal_generate",
            (uint16_t       )256,//任务堆栈大小
            (void*          )NULL,  
            (UBaseType_t    )1, //任务优先级
            (TaskHandle_t*  )&SignalGenerate_Handler);
vTaskStartScheduler();




int32_t sin_signal_generate()
{
        while(time * 1000 > tt && i < 63)
        {
                i++;
        }
        if(i == 64)
        {
                sin_signal = 0;
                time = 0;
                sin_signal_switch = 0;
        }
       
        sin_signal = sin(2 * pi * f * time);//正弦信号
        //sin_signal = sin(2 * pi * f * time) > 0?0:1;//阶跃信号
        time += 0.002;
       
        return E_OK;
}


然后就是使用sin_signal驱动电机，注意这里产生的信号是-1到1之间的，在驱动前应变换到合适的范围。编译好文件后，Keil进入debug模式，打开J-link并设置好目标——PID速度环输出和电机位置输出。注意J-Scope只能识别int整形数据，所以需要在程序里乘以1000然后用J-Scope读取并导入Matlab。


使用systemIdentification工具，选择 Import data > Time domain data ，在 input 和 output 选项中写上工作区的输入输出变量名，Starting time 和 Sample time 分别输入 0 和 0.002（这里的Sample time需要和J-Scope设置的一样）。点击 Import 将数据加入System Identification工具。


选择刚刚导入的数据，Estimate > State Space Model，得到系统传递函数如下：


ss1 =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)

  A =
              x1         x2         x3         x4
   x1     0.9951  -0.007623  0.0001774  -0.008001
   x2    0.05244     -0.267    -0.5505     0.6961
   x3  -0.005999     0.4216    -0.7989    -0.4194
   x4    0.04963     0.6398     0.2278      0.156

  B =
               u1
   x1   1.229e-09
   x2   2.354e-07
   x3  -1.292e-07
   x4   8.957e-08

  C =
               x1          x2          x3          x4
   y1  -8.353e+04       -1378      -372.7       207.4

  D =
       u1
   y1   0

  K =
               y1
   x1  -4.332e-06
   x2    0.000132
   x3   -7.29e-05
   x4  -5.293e-05

Name: ss1
Sample time: 0.002 seconds
  
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: estimate
   Number of free coefficients: 28
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.


Status:                                          
Estimated using N4SID on time domain data "cs".  
Fit to estimation data: 97.95% (prediction focus)
FPE: 30.51, MSE: 30.5






其中yaw_step是一个简单的阶跃信号，ss1是刚刚便是出来的方程。将单片机程序中的PID参数导入模型，得到模拟结果。
下图是我便是得到的结果及PID参数优化后得到的结果




可以看出，辨识出来的模型能较好的反应真实系统做出的响应，可以用于PID调参及后续射击系统建模。