<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE exportData>
<exportData ApplicationVersion="2.4.279-b11" ApplicationName="zWorkbench">
  <LibraryRefs/>
  <PrjItem InstanceFactoryInfo="ztools.FactoryTag.StPou" ID="1269" Name="АНР">
    <Property text="(*АНР-2: Автоматическая настройка ПИД по методу релейного эксперимента&#10;         с измерением двух периодов автоколебаний.&#10;&#10;  ПРИНЦИП РАБОТЫ:&#10;  ───────────────&#10;  1. В контур управления вместо ПИД вводится релейный элемент с гистерезисом H.&#10;     Релей переключает выход между (μ_пред + ΔMu) и (μ_пред - ΔMu).&#10;  2. Процесс входит в режим устойчивых автоколебаний. Измеряем 2 периода.&#10;  3. По форме колебаний вычисляем:&#10;       - Амплитуду выходной переменной  Ay&#10;       - Амплитуду управляющего сигнала Aμ (с учётом динамики ИМ)&#10;       - Модуль КЧХ объекта  Rob = Ay / Aμ&#10;       - Фазу  КЧХ объекта  Фob ≈ -π (условие автоколебаний)&#10;  4. Аппроксимируем объект моделью 1-го порядка с запаздыванием:&#10;       W(s) = Kob / (T1·s + 1) · e^(-β·T1·s)&#10;     Параметры Kob, T1, β находим из Rob, Фob, n=T1/β·T1 методом итераций.&#10;  5. Рассчитываем оптимальные Kp, Ti, Td по критерию минимума интегрального&#10;     отклонения с запасами устойчивости (ρs_op = 1.1, Фs_op = -70°).&#10;&#10;  СОСТОЯНИЯ (eState):&#10;  ────────────────────&#10;  0 — IDLE        Ожидание команды bStart или bpodstroika&#10;  1 — INIT        Инициализация переменных&#10;  2 — RELAY       Релейный режим + накопление измерений&#10;  3 — CALCULATE   Расчёт параметров модели и ПИД&#10;  4 — DONE        Готово: rKp, rTi, rTd действительны&#10;  9 — ERROR       Таймаут или ошибка (→ bReset для сброса)&#10;&#10;  ВХОДЫ:&#10;  ───────&#10;  bStart        — запуск полной автонастройки (ждём 2 периода)&#10;  bpodstroika   — режим подстройки (запуск без флага bBusy)&#10;  bEnable       — разрешение работы блока&#10;  bReset        — сброс в IDLE&#10;  rSetpoint     — уставка процесса&#10;  rProcessVar   — текущее значение процесса (PV)&#10;  rMuFeedback   — текущее значение управляющего сигнала (обратная связь)&#10;  rY0           — начальное (исходное) значение PV до пуска, дефолт 10.0&#10;  rTm           — время хода исполнительного механизма (ИМ), сек, дефолт 30.0&#10;  rMuMin/rMuMax — диапазон управляющего сигнала, дефолт 0..100 %&#10;  rDMu          — амплитуда релейного сигнала (полная), дефолт 30 %&#10;  rCycleTime    — время цикла вызова программы, сек (дефолт 0.1)&#10;  rKf           — коэффициент фильтра при расчёте Ti (дефолт 8.0)&#10;══════════════════════════════════════════════════════════════════════════════ *)&#10;PROGRAM ANR_2&#10;VAR_INPUT&#10;    bStart        : BOOL;               // Старт автонастройки&#10;    bReset        : BOOL;               // Сброс → IDLE&#10;    bEnable       : BOOL := TRUE;       // Разрешение работы блока&#10;    bpodstroika   : BOOL := FALSE;      // Режим подстройки (без bBusy)&#10;    rSetpoint     : REAL;               // Уставка процесса&#10;    rProcessVar   : REAL;               // Текущее значение PV&#10;    rMuFeedback   : REAL;               // Текущий управляющий сигнал (обратная связь ИМ)&#10;    rY0           : REAL := 10.0;       // Начальное значение PV (до автонастройки)&#10;    rTm           : REAL := 30.0;       // Время хода ИМ, сек&#10;    rMuMin        : REAL := 0.0;        // Минимум управляющего сигнала&#10;    rMuMax        : REAL := 100.0;      // Максимум управляющего сигнала&#10;    rDMu          : REAL := 30.0;       // Полная амплитуда релейного сигнала (ΔMu)&#10;    rCycleTime    : REAL := 0.1;        // Время цикла, сек&#10;    rKf           : REAL := 8.0;        // Коэффициент при расчёте Ti_op&#10;END_VAR&#10;&#10;VAR_OUTPUT&#10;    rControlOut   : REAL;               // Выход: управляющий сигнал (в режиме RELAY = рeлe)&#10;    rKp1           : REAL;              // Рассчитанный коэффициент Kp ПИД&#10;    rTi1           : REAL;              // Рассчитанная постоянная интегрирования Ti, сек&#10;    rTd1           : REAL;              // Рассчитанная постоянная дифференцирования Td, сек&#10;    bReady1        : BOOL;              // TRUE = расчёт завершён, параметры готовы&#10;    bBusy         : BOOL;               // TRUE = идёт эксперимент&#10;    eStatus       : BYTE;               // Код статуса (совпадает с eState)&#10;    rU11           : REAL;              // Нижний порог переключения реле (SP - H)&#10;    rU21           : REAL;              // Верхний порог переключения реле (SP + H)&#10;    period           : INT;             // Период автоколебаний, сек        // Найденное отношение n = T1 / (β·T1)&#10;END_VAR&#10;&#10;VAR&#10;    eState        : BYTE ;              // Текущее состояние автомата&#10;    rKp           : REAL;               // Рассчитанный коэффициент Kp ПИД&#10;    rTi           : REAL;               // Рассчитанная постоянная интегрирования Ti, сек&#10;    rTd           : REAL;               // Рассчитанная постоянная дифференцирования Td, сек&#10;    // ── Параметры гистерезиса и переключения реле ─────────────────────────────&#10;    rH            : REAL;               // Полуширина гистерезиса реле (H ≈ 4% от |SP-Y0|)&#10;    rU1           : REAL;               // Нижний порог переключения реле (SP - H)&#10;    rU2           : REAL;               // Верхний порог переключения реле (SP + H)&#10;    rMuPrev       : REAL;               // Предыдущее значение управляющего сигнала&#10;    bFirstSwitch  : BOOL;               // TRUE = первое переключение (половинная амплитуда)&#10;    rMuCurrent    : REAL;               // Текущий выход реле&#10; rAy           : REAL;               // Амплитуда колебаний PV (измерено)&#10;    rAmu          : REAL;               // Амплитуда 1-й гармоники управляющего сигнала&#10;    // ── Измерение длительности полупериодов ───────────────────────────────────&#10;    tEdgeTimer    : TON;                // Вспомогательный таймер (не используется активно)&#10;    rTon          : REAL;               // Суммарное время пребывания PV выше U2 (полупериод +)&#10;    rToff         : REAL;               // Суммарное время пребывания PV ниже U1 (полупериод -)&#10;    iPeriodCount  : INT := 0;           // Счётчик пройденных периодов (0→1→2→3)&#10;    rCurrentTime  : REAL := 0.0;        // Накопленное время (для возможного расширения)&#10;    bPrevAboveU2  : BOOL;               // Предыдущее состояние: PV &gt; U2 (для детекта фронта)&#10;    bRelayOff : BOOL; &#10;    bRelayOn  : BOOL; &#10;    // ── Накопители для расчёта амплитуды по методу дисперсии ─────────────────&#10;    // Амплитуда Ay вычисляется из среднеквадратичного отклонения ошибки:&#10;    //   σ² = &lt;e²&gt; - &lt;e&gt;²  →  Ay ≈ sqrt(2·σ²)&#10;    rS0_acc       : REAL := 0.0;        // ∫ e(t) dt  / Tn   → среднее смещение (C0)&#10;    rSc_acc       : REAL := 0.0;        // ∫ e²(t) dt / Tn   → среднеквадратичная ошибка&#10;    rSmu_acc      : REAL := 0.0;        // ∫ μ(t) dt  / Tn   → среднее значение μ&#10;&#10;    // ── Промежуточные переменные расчётного блока ─────────────────────────────&#10;    rError        : REAL;               // Текущая ошибка e = SP - PV&#10;    rC0           : REAL;               // Среднее смещение PV от уставки&#10;    rSm           : REAL;               // Скорость хода ИМ (% / сек)&#10;    rTs           : REAL;               // Время нарастания реле до ΔMu, сек&#10;    rB1           : REAL;               // π·Ts/Tn — нормированная длит. фронта ИМ&#10;    rB2           : REAL;               // π·Ton/Tn — нормированная длит. полупериода +&#10;    rS1           : REAL;               // sinc(B1) — ослабление из-за конечного фронта ИМ&#10;    rS2           : REAL;               // sin(B2)  — вклад асимметрии полупериода&#10;    rUnderRoot    : REAL;               // Подкоренное выражение при расчёте Ay&#10;    bReady        : BOOL;&#10;    // ── Промежуточные для решения уравнений (поиск Ω, β) ─────────────────────&#10;    rCval         : REAL;               // (Kob/Rob)² — входная величина для квадратного ур-я&#10;    rZval         : REAL;               // n² = (N_Fixed)²&#10;    rDval         : REAL;               // Дискриминант квадратного уравнения&#10;    rQval         : REAL;               // Корень: ω² = QVAL&#10;    rSm_tmp       : REAL;               // (не используется как rSm, временная)&#10;    rTs_tmp       : REAL;               // Временная переменная&#10;rTn : REAL; &#10;    rRob          : REAL;               // Модуль КЧХ объекта на частоте автоколебаний&#10;    rFob          : REAL;               // Фаза КЧХ объекта, рад&#10;    &#10;    rKob_calc     : REAL;               // Статический коэффициент усиления объекта Kob&#10;    rT1_calc      : REAL;               // Постоянная времени T1 модели объекта, сек&#10;    rN_Fixed      : REAL;  &#10;    // ── Параметры модели объекта ──────────────────────────────────────────────&#10;    // Используется модель: W(s) = Kob / [(T1·s+1)(n·T1·s+1)]&#10;    // (апериодическое звено 2-го порядка, разложение через β)&#10;    rKob          : REAL := 1.0;        // Статический коэффициент усиления объекта&#10;    rT1           : REAL := 1.0;        // Постоянная времени T1 (большая), сек&#10;    rBeta         : REAL;               // β = T2/T1 ∈ [0.05..1.0] (отношение постоянных)&#10;    rOmega        : REAL;               // Нормированная частота автоколебаний ω = Ω·T1&#10;&#10;    // ── Переменные для оптимального регулятора ────────────────────────────────&#10;    rOmega_op     : REAL;               // Нормированная оптимальная частота ω_op&#10;    rKs_op        : REAL;               // Модуль АЧХ регулятора на частоте ω_op&#10;&#10;    // ── Промежуточные для расчёта КЧХ регулятора ─────────────────────────────&#10;    rYs           : REAL;               // Среднее значение PV с учётом смещения&#10;    rMus          : REAL;               // Среднее значение μ&#10;    rAr           : REAL;               // Вещественная часть КЧХ регулятора&#10;    rBr           : REAL;               // Мнимая часть КЧХ регулятора&#10;    rRr_op        : REAL;               // Модуль КЧХ регулятора на ω_op&#10;    rFr_op        : REAL;               // Фаза КЧХ регулятора на ω_op&#10;    rFob_op       : REAL;               // Требуемая фаза объекта на ω_op&#10;    rCval2        : REAL;               // Промежуточная: Yn_Ti / (2π)&#10;    rZval2        : REAL;               // Промежуточная: α·Kf / Cval2&#10;    rXval         : REAL;               // Zval2²&#10;    rYval         : REAL;               // (Xval+1)²·Kf&#10;&#10;    // ── Итерации Ньютона (поиск ω_op) ────────────────────────────────────────&#10;    iIter         : INT;                // Счётчик итераций метода Ньютона&#10;    rXn           : REAL;               // Текущее приближение ω_op&#10;    rGn           : REAL;               // Значение функции G(ω) на шаге итерации&#10;    rGGn          : REAL;               // Производная G'(ω) на шаге итерации&#10;&#10;    // ── Аппроксимационные переменные (таблицы рекомендаций) ──────────────────&#10;    rAlpha        : REAL;               // α = Td/Ti (отношение постоянных Д и И)&#10;    rTi_op        : REAL;               // Оптимальное Ti (промежуточный результат)&#10;    rYn_alpha     : REAL;               // Интерполяция для α по n&#10;    rYn_Ti        : REAL;               // Интерполяция для Ti по β и n&#10;    rZ_beta       : REAL;               // Коэффициент интерполяции по β&#10;    rY01          : REAL;               // Точка интерполяции Ti при β=0.1&#10;    rY08          : REAL;               // Точка интерполяции Ti при β=0.8&#10;&#10;    // ── Итерации поиска β (REPEAT..UNTIL) ────────────────────────────────────&#10;    iNiter        : INT;                // Счётчик итераций корректировки n&#10;&#10;    // ── Таймаут безопасности ──────────────────────────────────────────────────&#10;    tTimeout      : TON;                // Таймер безопасности (блок автоколебаний)&#10;    tSafetyLimit  : TIME := T#800S;     // Максимальное время эксперимента, сек&#10;&#10;    // ── Математические константы ──────────────────────────────────────────────&#10;    PI_2          : REAL := 6.28318530717958;   // 2π&#10;    PI            : REAL := 3.14159265358979;   // π&#10;&#10;    // ── Коэффициенты аппроксимационных таблиц оптимального ПИД ───────────────&#10;    // Формулы вида: y = b + c/(x + a)  или  y = b - c/(x + a)&#10;    // Таблица α (отношение Td/Ti):&#10;    a03:REAL:=2.3107;  b03:REAL:=0.08498; c03:REAL:=0.7814;   // α_n  = b03 + c03/(n+a03)&#10;    a6 :REAL:=0.06518; b6 :REAL:=0.1425;  c6 :REAL:=0.01296;  // γ6  (вспомог. для α)&#10;    a40:REAL:=0.11897; b40:REAL:=0.03535; c40:REAL:=0.02705;  // γ40 (вспомог. для α)&#10;    d1 :REAL:=13.24;   d2 :REAL:=1.369;   d3 :REAL:=2.369;    // масштабы таблицы α&#10;    // Таблица Ti:&#10;    a01:REAL:=0.5833;  b01:REAL:=2.0133;  c01:REAL:=0.8444;   // Ti при β→0.1&#10;    a08:REAL:=0.6570;  b08:REAL:=7.052;   c08:REAL:=4.328;    // Ti при β→0.8&#10;    d4 :REAL:=0.1429;  d5 :REAL:=1.429;                        // масштабы по β&#10;    // Требования к устойчивости (точки оптимума на ГАС):&#10;    rRs_op    : REAL := 1.1;    // Желаемый модуль АЧХ разомкнутой системы&#10;    rGs_op    : REAL := -70.0;  // Желаемый запас по фазе, градусы (не используется напрямую)&#10;    rRrs_op   : REAL := 0.91;   // Целевой модуль КЧХ регулятора на ω_op&#10;    rFrs_op   : REAL := -2.24;  // Целевая фаза КЧХ регулятора на ω_op, рад&#10;END_VAR&#10;    rCurrentTime := rCurrentTime + rCycleTime;&#10;//  НАКОПЛЕНИЕ ВРЕМЕНИ (каждый скан при bEnable)&#10;IF NOT bEnable THEN&#10;    rCurrentTime := 0;&#10;    eState:=0;&#10;END_IF;&#10;//  ТАЙМЕР БЕЗОПАСНОСТИ: если эксперимент длится более tSafetyLimit → eState=9&#10;tTimeout(IN := bBusy, PT := tSafetyLimit);&#10;&#10;&#10;//  ГЛАВНЫЙ АВТОМАТ СОСТОЯНИЙ&#10;CASE eState OF&#10;&#10;0: (* IDLE — ожидание команды *)&#10;&#10;    bReady      := FALSE;&#10;    bBusy       := FALSE;&#10;    eStatus     := 0;&#10;    rControlOut := rMuFeedback;  // В режиме ожидания — прозрачная передача сигнала&#10;&#10;    IF bStart AND bEnable AND NOT tTimeout.Q THEN&#10;        // Полная автонастройка: выставляем флаг занятости&#10;        eState  := 1;&#10;        bBusy   := TRUE;&#10;        eStatus := 1;&#10;    ELSIF bpodstroika AND bEnable THEN&#10;        // Режим подстройки: без флага bBusy&#10;        eState  := 1;&#10;        eStatus := 1;&#10;    END_IF;&#10;&#10;1: (* INIT — инициализация всех переменных перед экспериментом *)&#10;&#10;    // Вычисляем гистерезис реле: H = 4% от расстояния |SP - Y0|, минимум 0.5&#10;    // Гистерезис нужен для помехозащищённости и формирования чётких переключений&#10;    rH  := 0.04 * ABS(rSetpoint - rY0);&#10;    IF rH &lt; 0.5 THEN rH := 0.5; END_IF;&#10;    rU1 := rSetpoint - rH;   // Нижний порог (при снижении PV ниже U1 → реле +ΔMu)&#10;    rU2 := rSetpoint + rH;   // Верхний порог (при росте PV выше U2 → реле -ΔMu)&#10;&#10;    // Инициализация состояния реле&#10;    rMuPrev      := rMuFeedback;   // Стартуем с текущего значения μ&#10;    rMuCurrent   := rMuFeedback;&#10;    bFirstSwitch := TRUE;          // Первое переключение — с половинной амплитудой&#10;&#10;    // Обнуляем счётчики периодов и накопители&#10;    iPeriodCount := 0;&#10;    rTon         := 0.0;&#10;    rToff        := 0.0;&#10;    rS0_acc      := 0.0;&#10;    rSc_acc      := 0.0;&#10;    rSmu_acc     := 0.0;&#10;    bRelayOff := FALSE;&#10;    bRelayOn  := FALSE;&#10;&#10;    rN_Fixed     := 10.0;    // Начальное приближение n (отношение постоянных времени)&#10;    bPrevAboveU2 := FALSE;   // PV ещё не пересекал U2&#10;&#10;    tEdgeTimer(IN := FALSE); // Сброс вспомогательного таймера&#10;&#10;    eState  := 2;&#10;    eStatus := 2;&#10;&#10;2: (* RELAY — релейный эксперимент + измерения *)&#10;//  Структура измерений:&#10;//  - Периоды считаются по нарастающим фронтам PV через U2 (снизу вверх).&#10;//  - Измерения Ton, Toff и интегралов ведутся ТОЛЬКО во 2-м периоде,&#10;//    т.к. 1-й период может быть несимметричным (переходный процесс).&#10;//  - По завершении 3-го периода (iPeriodCount &gt; 2) переходим в расчёт.&#10;&#10;    // Пересчитываем гистерезис каждый скан (уставка может меняться)&#10;    rH  := 0.04 * ABS(rSetpoint - rY0);&#10;    IF rH &lt; 0.5 THEN rH := 0.5; END_IF;&#10;    rU1 := rSetpoint - rH;&#10;    rU2 := rSetpoint + rH;&#10;&#10;    // ── Релейный элемент с гистерезисом ──────────────────────────────────────&#10;    IF bFirstSwitch THEN&#10;        // Первое переключение — вводим лишь половину амплитуды для мягкого старта&#10;        IF rProcessVar &lt; rSetpoint THEN&#10;            rMuCurrent := rMuPrev + 0.5 * rDMu;  // PV низко → греем/открываем&#10;        ELSE&#10;            rMuCurrent := rMuPrev - 0.5 * rDMu;  // PV высоко → охлаждаем/закрываем&#10;        END_IF;&#10;        bFirstSwitch := FALSE;&#10;    ELSE&#10;        // Нормальный режим: полная амплитуда ΔMu при выходе за гистерезис&#10;        IF rProcessVar &lt; rU1 THEN&#10;            rMuCurrent := rMuPrev + rDMu;   // PV ниже нижнего порога → +ΔMu&#10;        ELSIF rProcessVar &gt; rU2 THEN&#10;            rMuCurrent := rMuPrev - rDMu;   // PV выше верхнего порога → -ΔMu&#10;        END_IF;&#10;    END_IF;&#10;    // Ограничиваем выход диапазоном [MuMin..MuMax]&#10;    rMuCurrent := LIMIT(rMuMin, rMuCurrent, rMuMax);&#10;&#10;    // ── Счёт периодов: детект нарастающего фронта PV через U2 ────────────────&#10;    // Каждый раз, когда PV пересекает U2 снизу вверх — начался новый период&#10;    IF bRelayOff AND NOT bPrevAboveU2 THEN&#10;        iPeriodCount := iPeriodCount + 1;&#10;        bPrevAboveU2 := TRUE;&#10;    END_IF;&#10;    IF bRelayOn THEN&#10;        bPrevAboveU2 := FALSE;&#10;    END_IF;&#10;&#10;IF rProcessVar &gt; rU2 THEN&#10;    bRelayOff := TRUE;&#10;    bRelayOn  := FALSE;&#10;ELSIF rProcessVar &lt; rU1 THEN&#10;    bRelayOff := FALSE;&#10;    bRelayOn  := TRUE;&#10;END_IF;&#10;&#10;    // ── Измерение полупериодов Ton и Toff во 2-м периоде ─────────────────────&#10;    // Ton  = суммарное время, когда PV выше U2 (зона «включено»)&#10;    // Toff = суммарное время, когда PV ниже U1 (зона «выключено»)&#10;    // Накапливаем только во 2-м периоде, чтобы избежать переходного 1-го период&#10;    IF iPeriodCount = 2 THEN&#10;        IF bRelayOff THEN&#10;            rToff := rToff + rCycleTime;&#10;        END_IF;&#10;        IF bRelayOn THEN&#10;            rTon := rTon + rCycleTime;&#10;        END_IF;&#10;    END_IF;&#10;&#10;    // ── Накопление интегралов для расчёта амплитуды (2-й период) ─────────────&#10;    // ∫e dt  → rS0_acc  (для нахождения среднего смещения C0)&#10;    // ∫e² dt → rSc_acc  (для дисперсии → амплитуда Ay)&#10;    // ∫μ dt  → rSmu_acc (для оценки среднего μ → статический Kob)&#10;    IF iPeriodCount = 2 THEN&#10;        rError   := rSetpoint - rProcessVar;&#10;        rS0_acc  := rS0_acc  + rError            * rCycleTime;&#10;        rSc_acc  := rSc_acc  + rError * rError   * rCycleTime;&#10;        rSmu_acc := rSmu_acc + rMuFeedback        * rCycleTime;&#10;    END_IF;&#10;&#10;    // ── Завершение: 2 полных периода отработаны → переходим к расчёту ────────&#10;    IF iPeriodCount &gt; 2 THEN&#10;        rTn    := rTon + rToff;   // Полный период = Ton + Toff&#10;        eState := 3;&#10;        eStatus:= 3;&#10;    END_IF;&#10;&#10;    // ── Защита от зависания ───────────────────────────────────────────────────&#10;    IF tTimeout.Q THEN&#10;        eState  := 9;&#10;        eStatus := 5;&#10;    END_IF;&#10;&#10;    rMuPrev     := rMuCurrent;&#10;    rControlOut := rMuCurrent;&#10;&#10;3: (* CALCULATE — расчёт параметров объекта и ПИД *)&#10;//  Этапы расчёта:&#10;//  1.  Ay   — амплитуда PV из дисперсии накопленных интегралов&#10;//  2.  Aμ   — амплитуда 1-й гармоники управляющего сигнала (с учётом ИМ)&#10;//  3.  Rob  — модуль КЧХ объекта  Rob = Ay / Aμ&#10;//  4.  Фob  — фаза КЧХ объекта  Фob ≈ -π (условие Барклгаузена)&#10;//  5.  Kob  — статический коэффициент усиления объекта&#10;//  6-7. Ω, β — параметры модели объекта W(s) = Kob/[(T1s+1)(nT1s+1)]&#10;//  8.  T1   — постоянная времени объекта&#10;//  9.  α    — соотношение Td/Ti (аппроксимация по таблицам)&#10;//  10. Ti_op — оптимальная Ti (аппроксимация по таблицам)&#10;//  11. КЧХ регулятора на оптимальной частоте&#10;//  12. Kp, Ti, Td — итоговые параметры ПИД&#10;&#10;    // ── 1. Амплитуда Ay колебаний PV ─────────────────────────────────────────&#10;    // Из накопленных интегралов: C0 = &lt;e&gt;, σ² = &lt;e²&gt; - C0²&#10;    // Ay ≈ sqrt(2·σ²) — для синусоидального сигнала sqrt(2)·σ = амплитуда&#10;    IF rTn &gt; 0.0 THEN&#10;        rC0        := rS0_acc / rTn;                         // Среднее отклонение&#10;        rUnderRoot := 2.0 * (rSc_acc / rTn - rC0 * rC0);    // 2·σ²&#10;        IF rUnderRoot &gt; 0.0 THEN&#10;            rAy := SQRT(rUnderRoot);&#10;        ELSE&#10;            rAy := 0.001;   // Защита от нулевой амплитуды&#10;        END_IF;&#10;    ELSE&#10;        rAy := 0.001;&#10;    END_IF;&#10;    IF rAy &lt; 0.001 THEN rAy := 0.001; END_IF;&#10;&#10;    // ── 2. Амплитуда Aμ управляющего сигнала ─────────────────────────────────&#10;    // Учитывает конечное время нарастания ИМ (Ts = ΔMu / Sm):&#10;    //   Sm = 100 / Tm  — скорость хода ИМ (%/сек)&#10;    //   Aμ = (4·ΔMu/π) · sinc(π·Ts/Tn) · sin(π·Ton/Tn)&#10;    // sinc(B1) = sin(B1)/B1 — учёт сглаживания фронта ИМ&#10;    // sin(B2) — учёт несимметрии периода&#10;    rSm := 100.0 / rTm;                  // Скорость хода ИМ, %/сек&#10;    rTs := rDMu / rSm;                   // Время нарастания сигнала ИМ, сек&#10;    rB1 := (PI * rTs) / rTn;             // Нормированная длительность фронта&#10;    rB2 := (PI * rTon) / rTn;            // Нормированная длительность полупериода +&#10;    IF ABS(rB1) &gt; 0.001 THEN&#10;        rS1 := SIN(rB1) / rB1;           // sinc(B1): затухание из-за фронта ИМ&#10;    ELSE&#10;        rS1 := 1.0;                       // Мгновенный фронт: sinc → 1&#10;    END_IF;&#10;    rS2  := SIN(rB2);                     // sin(B2): асимметрия периода&#10;    rAmu := (4.0 * rDMu / PI) * rS1 * rS2;&#10;    IF rAmu &lt; 0.001 THEN rAmu := 0.001; END_IF;&#10;&#10;    // ── 3. Модуль КЧХ объекта ────────────────────────────────────────────────&#10;    // Rob = Ay / Aμ  (из условия автоколебаний: |W(jω)·N(Ay)| = 1,&#10;    //                 N(Ay) = Aμ/Ay — линеаризованный коэффициент реле)&#10;    rRob := rAy / rAmu;&#10;    IF rRob &lt; 0.001 THEN rRob := 0.001; END_IF;&#10;&#10;    // ── 4. Фаза КЧХ объекта ──────────────────────────────────────────────────&#10;    // Базовая фаза из условия автоколебаний: ∠W(jω) = -π&#10;    // Поправки: ASIN(2H/Ay) — из-за гистерезиса реле;&#10;    //           π·Ts/Tn      — из-за запаздывания ИМ;&#10;    //           π·((Ton/Tn)-0.5)² — из-за асимметрии периода&#10;    rFob := -PI;&#10;    IF rAy &gt; (2.0 * rH) THEN&#10;        rFob := rFob - ASIN(2.0 * rH / rAy);   // Поправка гистерезиса&#10;    END_IF;&#10;    IF rTs &gt; 0.0 THEN&#10;        rFob := rFob - PI * rTs / rTn;          // Поправка фронта ИМ&#10;    END_IF;&#10;    rFob := rFob - PI_2 * ((rTon/rTn) - 0.5) * ((rTon/rTn) - 0.5);  // Асимметрия&#10;&#10;    // ── 5. Статический коэффициент усиления Kob ──────────────────────────────&#10;    // Kob = |ΔY_ст / Δμ_ст| — из статической характеристики объекта&#10;    // Ys = SP + C0 — среднее значение PV (с учётом смещения от уставки)&#10;    // Mus = &lt;μ&gt;   — среднее управляющее воздействие за период&#10;    rYs  := rSetpoint + rC0;&#10;    rMus := rSmu_acc / rTn;&#10;    IF ABS(rMus) &gt; 0.001 THEN&#10;        rKob := ABS(rYs - rY0) / ABS(rMus);&#10;    ELSE&#10;        rKob := 1.0;   // Нет статики — единичный коэффициент&#10;    END_IF;&#10;    IF rKob &lt; 0.001 THEN rKob := 0.001; END_IF;&#10;    rKob_calc := rKob;&#10;&#10;    // ── 6-7. Поиск нормированной частоты Ω и коэффициента β ─────────────────&#10;    // Уравнение КЧХ модели на частоте ω: |W(jω)|·|G_relay| = 1&#10;    // Разворачивается в: (Kob/Rob)² = (1+ω²)·(1+(nω)²)&#10;    // Квадратное уравнение относительно x=ω²:&#10;    //   n²·x² + (1+n²)·x + (1 - (Kob/Rob)²) = 0&#10;    // β находим из фазового условия:&#10;    //   Фob = -ATan(ω) - ATan(n·ω) - β·ω    →  β = (-Фob - ATan(ω) - ATan(nω)) / ω&#10;    // Итерируем n, пока β ∈ [0.05..1.0]&#10;    rCval  := (rKob / rRob) * (rKob / rRob);  // (Kob/Rob)²&#10;    iNiter := 0;&#10;    REPEAT&#10;        rZval := rN_Fixed * rN_Fixed;           // n²&#10;        rDval := (1.0 + rZval) * (1.0 + rZval) - 4.0 * rZval * (1.0 - rCval);&#10;        IF rDval &gt;= 0.0 THEN&#10;            rQval := (SQRT(rDval) - (1.0 + rZval)) / (2.0 * rZval);&#10;        ELSE&#10;            rQval := 0.01;&#10;        END_IF;&#10;        rOmega := SQRT(ABS(rQval));            // ω = √(Q)&#10;&#10;        // β из фазового условия&#10;        IF rOmega &gt; 0.001 THEN&#10;            rBeta := (-rFob - ATAN(rOmega) - ATAN(rOmega * rN_Fixed)) / rOmega;&#10;        ELSE&#10;            rBeta := 0.5;&#10;        END_IF;&#10;&#10;        // Корректируем n, чтобы β попал в допустимый диапазон&#10;        IF rBeta &lt; 0.05 THEN&#10;            rN_Fixed := rN_Fixed + 1.0;   // β слишком мало → увеличиваем n&#10;        ELSIF rBeta &gt; 1.0 THEN&#10;            rN_Fixed := rN_Fixed - 1.0;   // β слишком велико → уменьшаем n&#10;        END_IF;&#10;        iNiter := iNiter + 1;&#10;        IF iNiter &gt; 100 THEN&#10;            EXIT;   // Защита от зависания&#10;        END_IF;&#10;    UNTIL (rBeta &gt;= 0.05) AND (rBeta &lt;= 1.0)&#10;    END_REPEAT;&#10;&#10;    // ── 8. Постоянная времени T1 ──────────────────────────────────────────────&#10;    // ω = Ω·T1 → T1 = ω·(Tn / 2π)&#10;    rT1 := rOmega * (rTn / PI_2);&#10;    IF rT1 &lt; 0.1 THEN rT1 := 0.1; END_IF;&#10;    rT1_calc := rT1;&#10;&#10;    // ── 9. Расчёт α = Td/Ti (аппроксимационные формулы) ──────────────────────&#10;    // α зависит от n и β. Используем два вспомогательных полинома γ6 и γ40,&#10;    // интерполированных по n, и взвешиваем по β.&#10;    rYn_alpha := b03 + c03 / (rN_Fixed + a03);        // Нормировочная функция по n&#10;    rAlpha    := (b40 + (c40 / (rBeta + a40))) * (d3 - rYn_alpha * d1)&#10;               - (b6  + (c6  / (rBeta + a6 ))) * (d2 - rYn_alpha * d1);&#10;&#10;    // ── 10. Расчёт Yn_Ti (оптимальное Ti через аппроксимацию) ────────────────&#10;    // Ti_opt зависит от β и n. Интерполируем между крайними значениями β=0.1 и β=0.8.&#10;    rZ_beta := d4 + d5 * (rBeta - 0.2);               // Весовой коэффициент по β&#10;    rY01    := b01 - c01 / ((1.0 / rN_Fixed) + a01);  // Ti при β ≈ 0.1&#10;    rY08    := b08 - c08 / ((1.0 / rN_Fixed) + a08);  // Ti при β ≈ 0.8&#10;    rYn_Ti  := rY01 + rZ_beta * (rY08 - rY01);        // Интерполяция&#10;&#10;    // ── 11. КЧХ оптимального регулятора на ω_op ──────────────────────────────&#10;    // Оптимальный ПИД-регулятор проектируется так, чтобы АЧХ разомкнутой системы&#10;    // на частоте ω_op обеспечивала заданные запасы устойчивости (rRrs_op, rFrs_op).&#10;&#10;    // Параметры КЧХ регулятора: C(jω) = 1 + 1/(jω·Ti) + jω·Td&#10;    //   → с учётом Kf (фильтр производной): C(jω) = (1 + Z·Kf)/(1 + Z)·Z - jCval2&#10;    rCval2 := rYn_Ti / PI_2;&#10;    rZval2 := rAlpha * rKf / rCval2;&#10;    rXval  := rZval2 * rZval2;&#10;    rYval  := (rXval + 1.0) * (rXval + 1.0) * rKf;&#10;    rAr    := 1.0 + (2.0 * rXval / rYval);&#10;    rBr    := (1.0 - rXval) * (rZval2 / rYval) - rCval2;&#10;    rRr_op := SQRT(rAr * rAr + rBr * rBr);   // Модуль КЧХ регулятора&#10;    rFr_op := ATAN(rBr / rAr);               // Фаза КЧХ регулятора, рад&#10;&#10;    // Требуемая фаза объекта на ω_op: Фob_op = Фrs_op - Fr_op&#10;    rFob_op := rFrs_op - rFr_op;&#10;&#10;    // ── Поиск ω_op методом Ньютона ───────────────────────────────────────────&#10;    // Решаем уравнение G(ω) = 0:&#10;    //   G(ω) = β·ω + ATan(ω) + ATan(n·ω) + Фob_op = 0&#10;    // G'(ω) = β + 1/(ω²+1) + n/(n²ω²+1)&#10;    rXn   := 0.2;     // Начальное приближение ω_op&#10;    iIter := 0;&#10;    REPEAT&#10;        rGn  := rBeta * rXn + ATAN(rXn) + ATAN(rXn * rN_Fixed) + rFob_op;&#10;        rGGn := rBeta&#10;              + 1.0 / (rXn * rXn + 1.0)&#10;              + rN_Fixed / (rN_Fixed * rN_Fixed * rXn * rXn + 1.0);&#10;        IF ABS(rGGn) &gt; 1.0E-6 THEN&#10;            rXn := rXn - rGn / rGGn;   // Шаг Ньютона&#10;        ELSE&#10;            EXIT;   // Производная близка к нулю — выходим&#10;        END_IF;&#10;        iIter := iIter + 1;&#10;        IF iIter &gt; 50 THEN&#10;            EXIT;   // Ограничение числа итераций&#10;        END_IF;&#10;    UNTIL ABS(rGn / rGGn) &lt; 0.01&#10;    END_REPEAT;&#10;    rOmega_op := ABS(rXn);&#10;    IF rOmega_op &lt; 0.001 THEN rOmega_op := 0.001; END_IF;&#10;&#10;    // Коэффициент усиления регулятора на ω_op:&#10;    //   Ks_op = Rrs_op / (Rr_op · |W(jω_op)|)&#10;    rKs_op := rRrs_op / (rRr_op *&#10;              SQRT((rOmega_op * rOmega_op + 1.0) *&#10;                   (rOmega_op * rOmega_op * rN_Fixed * rN_Fixed + 1.0)));&#10;&#10;    // ── 12. Параметры ПИД ────────────────────────────────────────────────────&#10;    // Kp = Ks_op / Kob  (коэффициент усиления регулятора / статика объекта)&#10;    // Ti = T1 · (2π/ω_op) / Yn_Ti&#10;    // Td = Ti · α&#10;    rKp := rKs_op / rKob;&#10;    rTi := (rT1 * (PI_2 / rOmega_op)) / rYn_Ti;&#10;    rTd := rTi * rAlpha;&#10;&#10;    // Ограничения — защита от физически бессмысленных значений&#10;    IF rKp &lt; 0.01   THEN rKp := 0.01;   END_IF;&#10;    IF rKp &gt; 1000.0 THEN rKp := 1000.0; END_IF;&#10;    IF rTi &lt; 0.1    THEN rTi := 0.1;    END_IF;&#10;    IF rTd &lt; 0.0    THEN rTd := 0.0;    END_IF;&#10;&#10;    eState  := 4;&#10;    eStatus := 4;&#10;&#10;4: (* DONE — параметры рассчитаны, ожидаем bReset *)&#10;    bReady  := TRUE;&#10;    bBusy   := FALSE;&#10;    eStatus := 4;&#10;    IF bReset THEN&#10;        bReady := FALSE;&#10;        eState := 0;&#10;    END_IF;&#10;&#10;&#10;9: (* ERROR — таймаут или ошибка, ожидаем bReset *)&#10;    bBusy   := FALSE;&#10;    eStatus := 5;&#10;    IF bReset THEN&#10;        eState := 0;&#10;    END_IF;&#10;&#10;END_CASE;&#10;&#10;//  ВЫХОДНОЙ СИГНАЛ&#10;IF eState = 2 THEN&#10;    rControlOut := rMuCurrent;       // Выход реле — идёт на объект&#10;ELSIF eState = 9 THEN&#10;    rControlOut := rMuMin;           // Авария — безопасный минимум&#10;ELSIF eState &gt;= 4 THEN&#10;    rControlOut := 0.0;              // После расчёта — нейтральный выход&#10;END_IF;&#10;&#10;//  СБРОС НАКОПИТЕЛЕЙ при bReset&#10;IF bReset THEN&#10;    rS0_acc      := 0.0;&#10;    rSc_acc      := 0.0;&#10;    rSmu_acc     := 0.0;&#10;    iPeriodCount := 0;&#10;    tEdgeTimer(IN := FALSE);&#10;END_IF;&#10;period:=iPeriodCount;&#10;rU11           := rU1;               // Нижний порог переключения реле (SP - H)&#10;rU21           := rU2;&#10;rKp1           := rKp;              &#10;rTi1           :=rTi;              &#10;rTd1:= rTd;&#10;bReady1 := bReady;               &#10;END_PROGRAM"/>
  </PrjItem>
</exportData>
