Модельна орієнтоване проектування на коліні, ідентифікація систем в MATLAB / Simulink

    
 
Привіт, Хабр!
 
Сьогодні я хочу показати простий приклад ідентифікації системи, грунтуючись на спостереженнях і експериментальних даних. Це перша і вкрай важлива ступінь у розробці системи управління пристроєм, який описати аналітично або неможливо, або дуже складно, або неохота . Для початку розглянемо метод «чорного ящика з котом », «сірий» і «білий» методи залишимо на наступний раз.
Цікавляться прошу під кат.
 
 

Устаткування і програмне забезпечення в наявності

 
     
  • Watchdog Watchduck:
     
    Виконує ту ж функцію, що й Watchdog timer, але стосовно до інженера — допомагає вийти з зависання.
     
     
  •  
  • Машина реального часу Speedgoat . Трохи простіше, ніж на заголовній ілюстрації, приблизно така:
     
     
    У даному ПК (а це по суті і є IBM сумісний ПК) присутні плати вводу-виводу з аналоговими, дискретними, та іншими входами і виходами, які я використовую як для збору даних, так і для управління ходом експерименту. Але найголовніше — я можу одним натисканням запускати на ній моделі Simulink в режимі жорсткого реального часу. Про такі машинках (хай вибачать мені зменшувально-пестливе подальше до Неї звернення :-) я багато ще можу розповісти. Якщо цікаво — можна виділити для цього окремий топік.
     
     
  •  
  • MATLAB / Simulink
     
    Інструментам я не зраджую. Як і в попередніх моїх постах, я буду використовувати Simulink для створення моделей, генерувати з них код і запускати його на цікавих залозках.
     
     
  •  
  • Двигун постійного струму з енкодером і зубчастої передачею від Pololu. Ось такий:
     
     
  •  
  • Arduino + Motor Shield від seeedstudio
     
    Для роботи з Arduino -під Simulink я використовую бібліотеку Embedded Coder Target for Arduino .
    Ця бібліотека дозволяє використовувати моделі Simulink для генерації коду повністю готового до компіляції під дану платформу. Можна використовувати також вбудовану підтримку Arduino в MATLAB / Simulink , але мені даний таргет більш до душі.
     
     
  •  
  • Макетна плата, провідники, харчування та інша дрібнота
     
  •  
Насправді, будь-який з представлених компонентів може бути замінений на відповідний функціональний аналог (крім качки, звичайно!). Але, так як робилося все «на коліні» — то й право вибору залишається за власником цього важливого суглоба.
 
У зібраному вигляді, якщо це можна назвати зібраністю, все виглядає так:
 
 
 

Постановка завдання

Є система — необхідно отримати її модель для подальшої роботи. Простіше нікуди.
Поведінка системи та її моделі повинно максимально збігатися в межах допустимих вхідних впливів. Метод «чорного ящика» має на увазі, що ми не будемо, або майже не будемо, приймати до уваги фізику процесів усередині системи, а будемо розглядати її з загальної точки зору теорії управління.
 
У деталі я вдаватися в цей раз не буду, спочатку все і так видно, на око.
Але слід врахувати, що за лаштунками залишаються архіважливі і надзвичайно цікава теми — планування експерименту (DOE ) і дисперсійний аналіз (ANOVA ), які заслуговують окремого розгляду.
 
Цільова система в деякому наближенні і з спрощенням складається з наступного набору підсистем:
 
     
  • Дискретна підсистема — Arduino
  •  
  • Безперервна механічна підсистема — двигун
  •  
  • Безперервна електрична підсистема — Motor Shield
  •  
 Вхід системи — задане значення ШІМ 0-100% .
Еесть можливість змінювати шпаруватість ШІМ сигналу, подаючи на вихід «аналогового» порту Arduino значення uint8 від 0 до 255. Необхідно дізнатися, як реагуватиме силова електроніка і двигун на цей вплив.
Тестовий вектор буде зашитий в контролер, складається він з 501 елемента і видається з дискретизацією 0.04с. Разом ~ 20 секунд тестового впливу.
Виглядає тестовий вектор так:
 
 
На виході системи я очікую побачити хоч щось подібне :-)
Для успішної ідентифікації тестовий вектор повинен покривати максимальну кількість режимів роботи системи, в яких вона буде, або може працювати. Метод «чорного ящика» не дозволяє повною мірою досліджувати аварійні режими зважаючи на обмеженість кількості експериментальних зразків. Я буду досліджувати систему в робочих режимах плавного пуску і старт-стоп.
 
 Вихід системи — швидкість обертання вала, кут повороту вала. З квадратурного енкодера, при обертанні вала двигуна, я отримую два меандра зі зміщенням по фазі на pi / 2. Їх обробляє апаратний декодер, який присутній на одній з конфігурованих плат вводу-виводу. При завантаженні моделі Simulink на машину реального часу відбувається «прошивка» ПЛІС (FPGA) для забезпечення необхідних входів / виходів. Тікі, з декодера множимо на коефіцієнт і перетворюємо в кут, а швидкість зміни їх кількості в швидкість. Енкодер двигуна робить 64 тріггірованія за один оборот, що дозволяє з хорошою точністю відслідковувати описані вище параметри. Основний алгоритм на машині реального часу встановлений мною на роботу з частотою 10 kHz — на виході я отримую експериментальні дані з періодом дискретизації 0,0001 с . Апаратна ж частина ПЛІС працює на частоті 33 MHz , що дозволяє, крім основного алгоритму декодування, встигати застосувати ще й дляантідребезговой алгоритми без шкоди точності вимірів на досить високих оборотах. Простими словами декодування квадратурного сигналу на FPGA описано в прекрасній статті fpga4fun -> QuadratureDecoder .
 
Дискретизація величин, тимчасові затримки, моменти інерції механічних деталей, тертя спокою, тертя ковзання, тертя кочення, механічний люфт, розбалансування, резонансні явища, комутація індуктивних ланцюгів, ковзний контакт, іскріння — не повний перелік факторів, з якими нам довелося б зіткнутися при спробі описати дану систему аналітично. А адже це — всього лише ардуінка з моторчиком !
Погоджуся, деякими нюансами цілком можна знехтувати, але на практиці, рішення знехтувати будь-яким фактором — не найпростіше рішення.
 
Так виглядає ШІМ з Arduino (жовтий) і напруга на чисто резистивної навантаженні, яка підключена до Шілд:
 
 
 
А ось графіки при підключеному двигуні в режимах, коли двигун ще не обертається, і вже обертається:
 
 
 
 
 
Комутація індуктивного навантаження з механічним контактом напівпровідниковими ключами — вкрай цікава тема для досліджень, але сьогодні ми уваги їй приділяти не будемо. Ми зріжемо по прямій — будемо розглядати тільки вхід і вихід системи, не вникаючи в фізику процесів.
 
 Моделі Simulink
На машині реального часу я запускаю ось таку модель:
 
 
 
Структура її вкрай проста. Є блоки для входів, виходів, загальні налаштування плат вводу-виводу, виведення інформації на дисплей машинки і проста внутрішня логіка перерахунку тиків від квадратурного енкодера в кут і швидкість.
 
Модель для отримання прошивки під Arduino виглядає наступним чином:
 
 
 
Таблиця істинності, яка відстежує знак сигналу і управляє входами Motor Shield
 
 
 
Код алгоритму можна подивитися під спойлером:
 C + +
//
// File: arduinoTestVector.cpp
//
// Code generated for Simulink model 'arduinoTestVector'.
//
// Model version                  : 1.229
// Simulink Coder version         : 8.6 (R2014a) 27-Dec-2013
// C/C++ source code generated on : Thu May 22 14:47:55 2014
//
// Target selection: arduino_ec.tlc
// Embedded hardware selection: Atmel->AVR
// Code generation objectives:
//    1. Execution efficiency
//    2. ROM efficiency
//    3. RAM efficiency
// Validation result: Not run
//
#include "arduinoTestVector.h"

// Constant parameters (auto storage)
const ConstParam_arduinoTestVector arduinoTestVector_ConstP = {
  // Expression: myData
  //  Referenced by: '<S1>/Constant1'

  { 0, 1, 3, 4, 5, 6, 8, 9, 10, 11, 13, 14, 15, 17, 18, 19, 20, 22, 23, 24, 26,
    27, 28, 29, 31, 32, 33, 34, 36, 37, 38, 40, 41, 42, 43, 45, 46, 47, 48, 50,
    51, 52, 54, 55, 56, 57, 59, 60, 61, 62, 0, -1, -3, -4, -5, -6, -8, -9, -10,
    -11, -13, -14, -15, -17, -18, -19, -20, -22, -23, -24, -26, -27, -28, -29,
    -31, -32, -33, -34, -36, -37, -38, -40, -41, -42, -43, -45, -46, -47, -48,
    -50, -51, -52, -54, -55, -56, -57, -59, -60, -61, -62, 0, 3, 5, 8, 10, 13,
    15, 18, 20, 23, 26, 28, 31, 33, 36, 38, 41, 43, 46, 48, 51, 54, 56, 59, 61,
    64, 66, 69, 71, 74, 77, 79, 82, 84, 87, 89, 92, 94, 97, 99, 102, 105, 107,
    110, 112, 115, 117, 120, 122, 125, 0, -3, -5, -8, -10, -13, -15, -18, -20,
    -23, -26, -28, -31, -33, -36, -38, -41, -43, -46, -48, -51, -54, -56, -59,
    -61, -64, -66, -69, -71, -74, -77, -79, -82, -84, -87, -89, -92, -94, -97,
    -99, -102, -105, -107, -110, -112, -115, -117, -120, -122, -125, 170, 170,
    170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170,
    170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170,
    170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170,
    170, 170, 170, -170, -170, -170, -170, -170, -170, -170, -170, -170, -170,
    -170, -170, -170, -170, -170, -170, -170, -170, -170, -170, -170, -170, -170,
    -170, -170, -170, -170, -170, -170, -170, -170, -170, -170, -170, -170, -170,
    -170, -170, -170, -170, -170, -170, -170, -170, -170, -170, -170, -170, -170,
    -170, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
    255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
    255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
    255, 255, 255, 255, 255, 255, -255, -255, -255, -255, -255, -255, -255, -255,
    -255, -255, -255, -255, -255, -255, -255, -255, -255, -255, -255, -255, -255,
    -255, -255, -255, -255, -255, -255, -255, -255, -255, -255, -255, -255, -255,
    -255, -255, -255, -255, -255, -255, -255, -255, -255, -255, -255, -255, -255,
    -255, -255, -255, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85,
    85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85,
    85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, -85, -85,
    -85, -85, -85, -85, -85, -85, -85, -85, -85, -85, -85, -85, -85, -85, -85,
    -85, -85, -85, -85, -85, -85, -85, -85, -85, -85, -85, -85, -85, -85, -85,
    -85, -85, -85, -85, -85, -85, -85, -85, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};

// Block signals (auto storage)
BlockIO_arduinoTestVector arduinoTestVector_B;

// Block states (auto storage)
D_Work_arduinoTestVector arduinoTestVector_DWork;

// Model step function
void arduinoTestVector_step(void)
{
  // local block i/o variables
  uint8_T rtb_MOTORSHIELD_IN1;
  uint8_T rtb_MOTORSHIELD_IN2;
  uint8_T rtb_Abs;
  boolean_T aVarTruthTableCondition_1;
  boolean_T aVarTruthTableCondition_2;
  boolean_T aVarTruthTableCondition_3;
  uint16_T rtb_Init;

  // Outputs for Enabled SubSystem: '<Root>/ENABLE_TEST' incorporates:
  //   EnablePort: '<S1>/Enable'

  // S-Function (sfunar_digitalInput): '<Root>/RUN_TEST'
  if (((uint8_T)digitalRead(((uint8_T)2U))) > 0) {
    if (!arduinoTestVector_DWork.ENABLE_TEST_MODE) {
      // InitializeConditions for UnitDelay: '<S1>/Unit Delay'
      arduinoTestVector_DWork.UnitDelay_DSTATE = false;

      // InitializeConditions for UnitDelay: '<S4>/FixPt Unit Delay2'
      arduinoTestVector_DWork.FixPtUnitDelay2_DSTATE = 1U;

      // InitializeConditions for UnitDelay: '<S4>/FixPt Unit Delay1'
      arduinoTestVector_DWork.FixPtUnitDelay1_DSTATE = 0U;
      arduinoTestVector_DWork.ENABLE_TEST_MODE = true;
    }

    // Switch: '<S4>/Init' incorporates:
    //   Constant: '<S4>/Initial Condition'
    //   Logic: '<S4>/FixPt Logical Operator'
    //   UnitDelay: '<S1>/Unit Delay'
    //   UnitDelay: '<S4>/FixPt Unit Delay1'
    //   UnitDelay: '<S4>/FixPt Unit Delay2'

    if (arduinoTestVector_DWork.UnitDelay_DSTATE ||
        (arduinoTestVector_DWork.FixPtUnitDelay2_DSTATE != 0)) {
      rtb_Init = 0U;
    } else {
      rtb_Init = arduinoTestVector_DWork.FixPtUnitDelay1_DSTATE;
    }

    // End of Switch: '<S4>/Init'

    // Selector: '<S1>/Selector' incorporates:
    //   Constant: '<S1>/Constant1'

    arduinoTestVector_B.Selector = arduinoTestVector_ConstP.Constant1_Value
      [(int16_T)rtb_Init];

    // Switch: '<S4>/Reset' incorporates:
    //   UnitDelay: '<S1>/Unit Delay'

    if (arduinoTestVector_DWork.UnitDelay_DSTATE) {
      // Update for UnitDelay: '<S4>/FixPt Unit Delay1' incorporates:
      //   Constant: '<S4>/Initial Condition'

      arduinoTestVector_DWork.FixPtUnitDelay1_DSTATE = 0U;
    } else {
      // Update for UnitDelay: '<S4>/FixPt Unit Delay1' incorporates:
      //   Constant: '<S1>/Constant'
      //   Sum: '<S1>/Add'

      arduinoTestVector_DWork.FixPtUnitDelay1_DSTATE = rtb_Init + 1U;
    }

    // End of Switch: '<S4>/Reset'

    // Update for UnitDelay: '<S1>/Unit Delay' incorporates:
    //   Constant: '<S3>/Constant'
    //   RelationalOperator: '<S3>/Compare'

    arduinoTestVector_DWork.UnitDelay_DSTATE = (rtb_Init == 500U);

    // Update for UnitDelay: '<S4>/FixPt Unit Delay2' incorporates:
    //   Constant: '<S4>/FixPt Constant'

    arduinoTestVector_DWork.FixPtUnitDelay2_DSTATE = 0U;
  } else {
    if (arduinoTestVector_DWork.ENABLE_TEST_MODE) {
      // Disable for Outport: '<S1>/TEST_OUT'
      arduinoTestVector_B.Selector = 0;
      arduinoTestVector_DWork.ENABLE_TEST_MODE = false;
    }
  }

  // End of S-Function (sfunar_digitalInput): '<Root>/RUN_TEST'
  // End of Outputs for SubSystem: '<Root>/ENABLE_TEST'

  // Truth Table: '<Root>/Truth Table'
  // Truth Table Function 'Truth Table': '<S2>:1'
  //  ClockWise
  // Condition '#1': '<S2>:1:10'
  aVarTruthTableCondition_1 = (arduinoTestVector_B.Selector > 0);

  //  CounterClockWise
  // Condition '#2': '<S2>:1:14'
  aVarTruthTableCondition_2 = (arduinoTestVector_B.Selector < 0);

  //  Stop
  // Condition '#3': '<S2>:1:18'
  aVarTruthTableCondition_3 = (arduinoTestVector_B.Selector == 0);
  if ((!aVarTruthTableCondition_1) && (!aVarTruthTableCondition_2) &&
      aVarTruthTableCondition_3) {
    // Decision 'D1': '<S2>:1:20'
    //  Stop
    // Action '3': '<S2>:1:48'
    rtb_MOTORSHIELD_IN1 = 0U;

    // Action '3': '<S2>:1:49'
    rtb_MOTORSHIELD_IN2 = 0U;
  } else {
    // Decision 'D2': '<S2>:1:20'
    if (aVarTruthTableCondition_1 && (!aVarTruthTableCondition_2) &&
        (!aVarTruthTableCondition_3)) {
      // Decision 'D2': '<S2>:1:22'
      //  ClockWise
      // Action '1': '<S2>:1:34'
      rtb_MOTORSHIELD_IN1 = 1U;

      // Action '1': '<S2>:1:35'
      rtb_MOTORSHIELD_IN2 = 0U;
    } else {
      // Decision 'D3': '<S2>:1:22'
      if ((!aVarTruthTableCondition_1) && aVarTruthTableCondition_2 &&
          (!aVarTruthTableCondition_3)) {
        // Decision 'D3': '<S2>:1:24'
        //  CounterClockWise
        // Action '2': '<S2>:1:41'
        rtb_MOTORSHIELD_IN1 = 0U;

        // Action '2': '<S2>:1:42'
        rtb_MOTORSHIELD_IN2 = 1U;
      } else {
        // Decision 'D4': '<S2>:1:24'
        // Decision 'D4': '<S2>:1:26'
        //  Default
        //  None
        // Action '4': '<S2>:1:55'
        rtb_MOTORSHIELD_IN1 = 0U;

        // Action '4': '<S2>:1:56'
        rtb_MOTORSHIELD_IN2 = 0U;
      }
    }
  }

  // End of Truth Table: '<Root>/Truth Table'

  // S-Function (sfunar_digitalOutput): '<Root>/MOTORSHIELD_IN1'
  digitalWrite(((uint8_T)8U), rtb_MOTORSHIELD_IN1);

  // S-Function (sfunar_digitalOutput): '<Root>/MOTORSHIELD_IN2'
  digitalWrite(((uint8_T)11U), rtb_MOTORSHIELD_IN2);

  // Abs: '<Root>/Abs'
  if (arduinoTestVector_B.Selector < 0) {
    rtb_Abs = (uint8_T)-arduinoTestVector_B.Selector;
  } else {
    rtb_Abs = (uint8_T)arduinoTestVector_B.Selector;
  }

  // End of Abs: '<Root>/Abs'

  // S-Function (sfunar_analogOutput): '<Root>/SPEEDPIN_A'
  analogWrite(((uint8_T)9U), rtb_Abs);
}

// Model initialize function
void arduinoTestVector_initialize(void)
{
  // Registration code

  // block I/O
  (void) memset(((void *) &arduinoTestVector_B), 0,
                sizeof(BlockIO_arduinoTestVector));

  // states (dwork)
  (void) memset((void *)&arduinoTestVector_DWork, 0,
                sizeof(D_Work_arduinoTestVector));

  // S-Function (sfunar_digitalInput): <Root>/RUN_TEST
  pinMode(((uint8_T)2U), INPUT);

  // InitializeConditions for Enabled SubSystem: '<Root>/ENABLE_TEST'
  // InitializeConditions for UnitDelay: '<S1>/Unit Delay'
  arduinoTestVector_DWork.UnitDelay_DSTATE = false;

  // InitializeConditions for UnitDelay: '<S4>/FixPt Unit Delay2'
  arduinoTestVector_DWork.FixPtUnitDelay2_DSTATE = 1U;

  // InitializeConditions for UnitDelay: '<S4>/FixPt Unit Delay1'
  arduinoTestVector_DWork.FixPtUnitDelay1_DSTATE = 0U;

  // End of InitializeConditions for SubSystem: '<Root>/ENABLE_TEST'

  // S-Function (sfunar_digitalOutput): <Root>/MOTORSHIELD_IN1
  pinMode(((uint8_T)8U), OUTPUT);

  // S-Function (sfunar_digitalOutput): <Root>/MOTORSHIELD_IN2
  pinMode(((uint8_T)11U), OUTPUT);

  // S-Function (sfunar_analogOutput): <Root>/SPEEDPIN_A
  pinMode(((uint8_T)9U), OUTPUT);
}

//
// File trailer for generated code.
//
// [EOF]
//
Весь проект можна скачати тут .
 
Модель, яку я буду використовувати для перевірки результатів:
 
 
 
Всі моделі можна скачати тут .
«Відкрив і запустив» може не спрацювати, є нюанси. З питаннями в личку, або в коментарі.
 
 Відео
  
 Результати
Ось, що вийшло в результаті:
 
 
 
Синій — це вхідний сигнал з межами від -255 до +255 . У перерахунку на вольти буде приблизно від -8,5 до 8,5 .
Зелений — швидкість обертання вала двигуна.
 
Бачимо тимчасову затримку, відсутність обертання, або вкрай мале обертання при подачі ШІМ менше 25% . Також спостерігаємо класичний апериодический перехідний процес.
 
А от інформація, яка виводиться під час проведення експерименту на дисплей, підключений до машини реального часу:
 
 
 
 Ідентифікація
Одна з підсистем, а саме двигун, може бути описаний наступною передавальної функцією:
 
Так нам говорить курс ТАУ, з неї ми і почнемо.
 
У складі MATLAB є прекрасний інструмент для ідентифікації систем — System Identification Toolbox . Він доступний як у вигляді графічного інтерфейсу, так і у вигляді набору функцій, які можна використовувати в скриптах MATLAB . Розглянемо спершу роботу в графічному інтерфейсі.
 
Імпортуємо дані, отримані в ході експерименту і тестове вплив:
 
 
Після імпорту нам стають доступні інструменти попередньої обробки даних. Ми маємо можливість, наприклад, розділити масиви даних на дві частини — для ідентифікації і окремо для верифікації. Але це для більш складних випадків, давайте ж швидше перейдемо до справи!
 
Вибираємо зі списку «Estimate» пункт Transfer Function Models .
 
Отримуємо передавальну функцію:
 
 
Давайте порівняємо поведінку отриманої функції й системи:
 
 
Синім відображений відгук передавальної функції на вихідний тестовий вектор. Видно, що дану функцію можна використовувати як модель вихідної системи тільки в режимі старт-стоп. І не дивно — нелінійне поведінка, яке присутнє в системі, не може бути описано таким способом.
 
Просуваємося далі, використовуємо шаблон моделі Хаммерштейна-Вінера , яка дає можливість описувати нелінійне поведінка системи:
 
 
 
В якості вхідної нелінійності виберемо Dead Zone — відсутність реакції системи на вхідний сигнал менш певного порогового значення. Такий тип нелінійності має взяти на себе опис тертя спокою і впливу постійних магнітів, яке має місце в системі.
Інші параметри за замовчуванням, тиснемо Estimate .
 
Порівнюємо отриманий результат з експериментальними даними:
 
 
Виглядає набагато краще! Не ідеально — помітно деяке аномальна поведінка, але ж ми тільки злегка торкнулися ідентифікації нелінійних систем.
 
Тепер нам залишається тільки експортувати дану модель в Simulink , де ми можемо приступити до розробки системи управління. Але про це в наступному топіку!
 
Дякую за увагу, сподіваюся було цікаво?
І ще раз дякую за відповіді на невелике опитування:
  

    
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.