Controle preditivo em quatro passos

A cada amostra \(k\) (a cada 1 ms):

1. Prever: com o modelo, calcular como a saída vai evoluir nos próximos instantes para um plano de movimentos de \(u\);
2. Avaliar: medir quão bom é o plano com uma função custo \(J\) — erro em relação à referência e esforço do atuador;
3. Escolher: ajustar o plano para deixar \(J\) o menor possível;
4. Aplicar só o primeiro movimento do plano — e, na amostra seguinte, repetir tudo com a nova medida.

O que o controlador “vê” em um instante

No instante \(k\) (linha vertical): o passado é medido; o futuro é previsto para o plano escolhido. O plano não precisa ser perfeito — só o 1º movimento é aplicado, e o plano é refeito na amostra seguinte.

O modelo ARX identificado

Modelo da identificação do Sallen-Key, \(T_s\) = 1 ms, identificado com a PRBS:

\[ y(k) = a_1\,y(k-1) + a_2\,y(k-2) + b_1\,u(k-1) + b_2\,u(k-2) + c \]

Prever as variações

Subtraindo a equação em \(k-1\) da equação em \(k\), o termo constante \(c\) some:

\[ \Delta y(k) = a_1\,\Delta y(k-1) + a_2\,\Delta y(k-2) + b_1\,\Delta u(k-1) + b_2\,\Delta u(k-2) \]

com \(\Delta y(k) = y(k) - y(k-1)\) e \(\Delta u(k) = u(k) - u(k-1)\).

• O controlador prevê quanto a saída vai mudar a partir da medida atual, e soma as mudanças a \(y(k)\);
• Não precisa do ponto de operação nem do termo \(c\);
• Os movimentos do plano já são variações de \(u\): \(\Delta u(k), \Delta u(k+1), \ldots\)

Por que as variações? Erro de modelo

A função custo

\[ J = \sum_{j=1}^{N_p}\big[r(k) - \hat y(k+j)\big]^2 \;+\; \lambda\sum_{i=0}^{N_u-1}\Delta u(k+i)^2 \]

• 1º termo: o erro previsto entre a referência e a saída, nos próximos \(N_p\) instantes;
• 2º termo: o tamanho dos movimentos de \(u\) — o esforço do atuador;
• \(\lambda\) pondera um contra o outro:
  • \(\lambda\) pequeno: o controlador corre atrás da referência com movimentos grandes;
  • \(\lambda\) grande: movimentos suaves, resposta mais lenta.

\(N_p\) e \(N_u\)

• \(N_p\) — horizonte de previsão: quantas amostras à frente o controlador olha;
• \(N_u\) — movimentos planejados: quantos valores futuros de \(\Delta u\) ele escolhe; depois deles, \(u\) fica constante até o fim do horizonte;
• O plano tem só \(N_u\) incógnitas: \(\Delta u(k), \ldots, \Delta u(k+N_u-1)\);
• Escolha: \(N_p\) = 60 (60 ms) e \(N_u\) = 4.

Como escolher \(N_p\)

• A planta leva ≈ 100 ms para acomodar sozinha (constante de tempo ≈ 26 ms);
• \(N_p\) deve cobrir a parte importante da resposta da malha fechada (≈ 50 ms): 60 ms basta;
• Horizonte muito curto: o controlador “não enxerga” o fim da subida e exagera.

Como escolher \(N_u\)

• Mais movimentos dão mais liberdade, mas o ganho satura rápido: 4 bastam;
• Cada incógnita a mais custa tempo em cada iteração do gradiente.

O peso \(\lambda\)

Escolha: \(\lambda\) = 3 — acomodação ≈ 50 ms sem sobressinal (a especificação do planejamento e dos ensaios P e PI).

O método do gradiente

• Com \(N_u\) = 4, \(J\) é uma função de 4 variáveis com a forma de uma tigela: tem um único ponto mais baixo;
• O gradiente \(\nabla J\) aponta para onde \(J\) sobe mais depressa; um passo no sentido contrário desce: \[ \Delta\mathbf u \;\leftarrow\; \Delta\mathbf u - \text{passo}\cdot\nabla J \]
• A derivada de \(J\) em relação a cada movimento usa a resposta ao degrau \(s_j\) do modelo (quanto \(\hat y(k+j)\) muda por unidade de \(\Delta u(k)\)), calculada uma vez na inicialização: \[ \frac{1}{2}\frac{\partial J}{\partial \Delta u(k+i)} = -\sum_{j} s_{j-i}\,\big[r - \hat y(k+j)\big] + \lambda\,\Delta u(k+i) \]

O gradiente passo a passo

Exemplo com 2 movimentos, no instante em que a referência sobe: cada iteração desce a tigela. Os primeiros passos fazem quase todo o trabalho; o resto é ajuste fino.

Quantas iterações por amostra?

Cada iteração simula o modelo 60 passos e calcula 4 derivadas: ≈ 490 multiplicações ≈ 2 000 ciclos ≈ 10 µs (≈ 4 ciclos por multiplicação, medido no SK08).

O método é viável?

• Convergir completamente nas trocas da referência levaria centenas de iterações: não cabe em 1 ms;
• Mas não é preciso convergir: com 50 iterações por amostra (≈ 0,5 ms) a resposta fica a poucos milissegundos da ótima — os primeiros passos fazem quase todo o trabalho;
• Cada amostra parte do plano da amostra anterior: a melhoria continua no instante seguinte;
• Número fixo de iterações → tempo de cálculo constante: bom para um laço de 1 ms;
• É o preço de resolver numericamente a cada amostra; a medida na placa confirma (ou não) a estimativa.

Ruído de medida

Parâmetros no main.c

const ParametrosMPC p = {
    .modelo = {           // ARX, Ts = 1 ms (identificacao do Sallen-Key)
        .a1 = 1.60918827f, .a2 = -0.62250084f,
        .b1 = 0.01096776f, .b2 = 0.00237673f,
    },
    .Np       = 60,       // horizonte de previsao: 60 ms
    .Nu       = 4,        // movimentos planejados
    .lambda   = 3.0f,     // peso de du
    .nIter    = 50,       // iteracoes do gradiente por amostra
    .uMin     = 0.0f,     // limites do atuador [V]
    .uMax     = 1.7f,
};

Na inicialização: a resposta ao degrau \(s_j\) do modelo e o tamanho do passo do gradiente.

Previsão em C

Simula o modelo nas variações, a partir das medidas, para o plano dU:

static void Preve(void)
{
    float dyA = y0m - y1m, dyB = y1m - y2m;   // dy(k), dy(k-1) medidos
    float duAnt = u1 - u2, yy = y0m;          // du(k-1); parte de y(k)
    for (j = 0U; j < Np; j++) {
        du = (j < Nu) ? dU[j] : 0.0f;         // depois de Nu, u constante
        dy = a1 * dyA + a2 * dyB + b1 * du + b2 * duAnt;
        duAnt = du;  dyB = dyA;  dyA = dy;
        yy += dy;  yp[j] = yy;                // soma as variacoes
    }
}

Uma iteração do gradiente em C

static void IteracaoGradiente(float r)
{
    Preve();                                    // y^ para o plano atual
    for (i = 0U; i < Nu; i++) {                 // derivada de J/2 em du_i
        g = lambda * dU[i];
        for (j = i; j < Np; j++)
            g -= s[j - i] * (r - yp[j]);        // s: resposta ao degrau
        grad[i] = g;
    }
    for (i = 0U; i < Nu; i++)
        dU[i] -= passo * grad[i];               // desce a tigela
}

Uma amostra do MPC em C

float MPC(float r, float y)
{
    y2m = y1m;  y1m = y0m;  y0m = y;            // medidas
    for (i = 0U; i + 1U < Nu; i++) dU[i] = dU[i + 1U];
    dU[Nu - 1U] = 0.0f;                         // parte do plano anterior
    for (it = 0U; it < nIter; it++) IteracaoGradiente(r);
    u = u1 + dU[0];                             // so o 1o movimento
    if (u < uMin) u = uMin;                     // saturacao do atuador
    if (u > uMax) u = uMax;
    u2 = u1;  u1 = u;
    return u;
}

Medição do tempo na placa

O CPU Timer 1 conta ciclos de CPU (5 ns); a struct desempenho aparece no Watch:

Por amostra, tV[k] grava o tempo de cálculo em µs.

Cuidado: Device_cal e -O2

• O projeto é compilado com otimização (-O2) para medir tempos realistas;
• O Device_init chama uma rotina de calibração da TI que altera um registrador da CPU (XAR1) sem restaurá-lo;
• Com -O2, o compilador guardava ali o endereço da rotina de amostragem: o programa travava na primeira amostra (Interrupt_illegalOperationHandler);
• Correção, em placa.c: Ts e o endereço da rotina ficam em variáveis volatile durante o Device_init.

O que gravar e observar

• Gravar o SK07, Resume, esperar ~2 s (1 s de acomodação + 0,8 s de ensaio);
• No gráfico do CCS, Capacity = 800; exportar rV, uV, yV e tV para dados/DadosMPCConvencional01.csv;
• Anotar da desempenho: calculoMax, lacoMax, ocupacaoMax e estouros;
• Comparar com a previsão: acomodação ≈ 53 ms, sem sobressinal; cálculo ≈ 0,5 ms (estimado);
• Variar nIter (10, 20, 50, 100) e ver o efeito na resposta e no tempo — até onde cabe em 1 ms?

Referências

• Camacho, E. F.; Bordons, C. Model Predictive Control. 2ª ed. Springer, 2007.
• Maciejowski, J. M. Predictive Control with Constraints. Prentice Hall, 2002.
• Rossiter, J. A. Model-Based Predictive Control: A Practical Approach. CRC Press, 2003.
• Boyd, S.; Vandenberghe, L. Convex Optimization. Cambridge University Press, 2004 (método do gradiente).

Firmware: ProjetosSK/SK07MPCNumerico/; réplica em Python: analise/convencional.py.