Saltar al contenido
Microcontroladores (PIC)

Predictor de Smith con PIC

Hola controleras y controleros, hoy vamos a aprender como implementar un Predictor de Smith en un microcontrolador PIC, utilizando como controlador primário un PID. Todo será implementado en el compilador CCS C.

Antes que nada te hago la invitación para 👉 Ver el Curso de PIC completo.

Contenidos mostrar

Predictor de Smith

El predictor de Smith es una estructura de control ampliamente utilizada en la industria para poder realizar control de procesos que presentan un retardo dominante. La estructura se presenta a continuación:

Predictor de Smith
5. Predictor de Smith

Donde P(s) es el proceso (Horno), C(S) será el controlador primário (En este caso usaremos un PID), Gn(s) es el modelo del proceso (función de transferencia del horno), Ln es el retardo del horno.

Mira la siguiente entrada para entender el funcionamiento del 👉 Predictor de Smith.

Predictor de Smith con PIC

Para la implementación de la estructura del predictor de Smith dentro del microcontrolador, tendremos que recurrir a la representación discreta de la estructura.

Predictor de Smith Discreto

Y para eso vamos a tomar la función de transferencia que obtuvimos de la Identificación del Horno de Temperatura.

G_n(s)=\dfrac{2.3553}{370s+1}e^{-30s}

Vamos a discretizar esa función de transferencia con un tiempo de muestreo de T=10s usando el retenedor de orden cero ZOH. Donde obtendremos una función de transferencia discreta de la siguiente forma:

G_n(z)=\dfrac{b_0z+b_1}{z+a_1}z^{-d}

G_n(z)=\dfrac{0.0628}{z-0.9733}z^{-3}

Implementando un controlador PID discreto, podemos proceder a montar nuestra estructura del predictor de Smith en un microcontrolador PIC.

Predictor de Smith en un Horno con PIC

Ecuaciones en diferencias

Para poder escribir la función de transferencia discreta dentro del microcontrolador PIC, se debe transformar dicha función en una ecuación en diferencias.

G_n(z)=\dfrac{\hat{y}(k)}{u(k)}=\dfrac{0.0628}{z-0.9733}z^{-3}

Se lleva a potencias negativas en z

\dfrac{\hat{y}(k)}{u(k)}=\dfrac{0.0628z^{-1}}{1-0.9733z^{-1}}z^{-3}

Se resuelve para obtener la salida del modelo

\hat{y}(k)(1-0.9733z^{-1})=0.0628z^{-1}u(k)z^{-3}

\hat{y}(k)-0.9733\hat{y}(k)z^{-1}=0.0628z^{-1}u(k)z^{-3}

\hat{y}(k)-0.9733\hat{y}(k-1)=0.0628u(k-1-3)

\hat{y}(k)=0.0628u(k-1-3)+0.9733\hat{y}(k-1)

Con eso ya podemos calcular el resultado de la salida del modelo (función de transferencia) del horno en lenguaje C.

//Salida del modelo nominal Pn(z)
yE= b0*u[kT-d]+b1*u[kT-1-d]-a1*yE_1;

Notemos que

Salida del modelo

Y que u es un vector. En este caso de 15 posiciones, donde la ultima posición corresponde al tiempo actual (presente) y las demás posiciones corresponden a muestras en el pasado.

Vector U con retardo

Implementación del Predictor de Smith en un PIC 18F4550

El controlador primário del PIC, sera el controlador PID explicado en la entrada Control de Temperatura de un Horno.

El circuito a ser implementado, será el mismo empleado en la entrada anterior de la identificación del horno.

Donde se ve que se usa 2 XBEE para realizar una comunicación inálambrica.

La respuesta del controlador sobre el horno

El código del PIC en CCS C y el código en MATLAB para la interfaz gráfica pueden ser descargados a continuación, solo basta con que compartas el contenido de este post con cualquiera de los siguientes botonos, de esa forma ayudas a que más personas aprendan sobre esto y contribuyes a mejorar la educación.

Otros Controladores implementados en PIC

Si deseas saber como implementar otro tipo de controladores como diversas configuraciones del PID, logica Fuzzy, Espacio de estados, RST, etc, te invito a que te matricules a mi curso de Control en Sistemas Embebidos, y puedes seleccionar el micro que más te guste (PIC o Arduino):

  • Curso de Sistemas de Control en Dispositivos Microcontrolados en UDEMY (PIC y ARDUINO)
  • Certificado de Aprobación una vez finalices el Curso
  • DESCUENTO si accedes directamente con los siguientes botones de acceso.
  • NOTA: Si buscas el curso directamente en UDEMY o si lo adquieres en otra plataforma distintas a las mostradas anteriormente NO OBTENDRÁS NINGUN DESCUENTO sobre el valor final del Curso.

Código en CCS C (PIC C) Predictor de Smith

#include <18F4550.h>
#device ADC=10

#FUSES NOWDT                    //No Watch Dog Timer
#FUSES WDT128                   //Watch Dog Timer uses 1:128 Postscale
#FUSES NOBROWNOUT               //No brownout reset
#FUSES NOLVP                    //No low voltage prgming, B3(PIC16) or B5(PIC18) used for I/O
#FUSES NOXINST                  //Extended set extension and Indexed Addressing mode disabled (Legacy mode)
#FUSES HS

#use delay(crystal=20000000)
#use rs232(baud=9600,parity=N,xmit=PIN_C6,rcv=PIN_C7,bits=8)

#define TC_CLK     PIN_B1  
#define TC_CS      PIN_C0 
#define TC_DATA    PIN_B0



#include "max6675.c"
#include <stdlib.h>
#include <string.h>

long deg=0;
int muest=0;
int16 control=0;
float R=0;//Referencia de 150 °C por defecto
float yM=0,e=0.0,e_1=0.0,e_2=0.0,T=10;
float kp,ti,td,q0,q1,q2;
float k=2.3553 ,tao=370 ,theta=30;
float TsMA,Wn,P1,P2;
float b0=0,b1=0.0628,a1=-0.9733; //Parametros de la Funcion de transferencia discreta
int d=3-1; //Retardo discreto (d-1)
float u[15]; //Ley de control (debe ser mayor que el retardo + el numero de ceros de la FT)
int kT=14; //Tiempo discreto actual (Ultima posición del vector u)
float yE=0,yE_1=0; //Salida del modelo interno
float ep=0; //Error de prediccion
float yR=0,yR_1=0; //Salida del modelo Rapido (sin retardo)
float yp=0; //Salida  predicha


//Define la interrupción por recepción Serial
#int_RDA
void RDA_isr()
{
 int i=0,ini=0,fin=0;
 char dat[5];
 char degC[5];
 
 // Almacena 10 datos leidos del serial
  for(i=0;i<5;i++){
     dat[i]=getc();
     putc(dat[i]);
    }



//Busco el valor del escalon en los datos recibidos
for(i=0;i<5;i++){
  if(dat[i]=='S'){
    ini=i+1;
    i=10;
  }
 }
 for(i=ini;i<5;i++){
  if(dat[i]=='$'){
    fin=i-1;
    i=10;
  }
 }
 if(ini!=0 && fin!=0){
    // salvo en degC el caracter con el escalon
    for(i=ini;i<=fin;i++){
     degC[i-ini]=dat[i];
    }
    
     deg = atol(degC); //Convierte el String en un valor numerico
     
     // valida que el escalón esté ente 0 y 300 grados celcius
     
     if(deg>1000)
         deg=1000;
     if(deg<0)
         deg=0;
    
       R=deg;
 }
}
//! 

/*===========================================================================*/
/*=======================    FUNCION DEL CONTROL PID  =======================*/
/*===========================================================================*/
void PID(void)
{
    //char msg[32];
    e=1*(R-yp);
    // Controle PID
      u[kT] = u[kT-1] + q0*e + q1*e_1 + q2*e_2; //Ley del controlador PID discreto
      
//!     sprintf(msg,"e=%3.2f u=%3.2f",e,u); 
//!      printf("%s",msg);
      
      
    if (u[kT] >= 1000.0)        //Saturo la accion de control 'uT' en un tope maximo y minimo
     u[kT] = 1000.0;
    
    if (u[kT] <= 0.0)
     u[kT] = 0.0;
     
     control=u[kT];
     
     //Actualizo los valores pasados
     e_2=e_1;
     e_1=e;
     
     //La accion calculada la transformo en PWM
     set_pwm1_duty(control);
     output_toggle(PIN_B2);  //LED verifica el tiempo de muestreo correcto
}


void main() 
{ 
   int i;
   char msg[32]; 
   delay_ms(50);      //allow oscillator to stabilise 
   setup_timer_2(t2_div_by_4,249,1);   //Configuracion de Timer 2 para establecer frec. PWM a 1kHz
   setup_ccp1(ccp_pwm);                //Configurar modulo CCP1 en modo PWM
   set_pwm1_duty(0);                   //Dejo en cero la salida PWM
   enable_interrupts(INT_RDA);
   enable_interrupts(GLOBAL);
   set_tris_b(0);
   
   
   //*************************************************************************//
   //*************  DISEÑO POR ASIGNACIÓN DE 2 POLOS REALES   ****************//
   //*************************************************************************//
   
   TsMA=600;                    //Tiempo deseado en Lazo Cerrado    
   Wn=4/(TsMA);               //Frecuencia natural del sistema
   
   //Ubicación de 2 Polos reales
   P1=2*Wn;
   P2=Wn*Wn;
   
   kp=(P1*tao-1)/k;        //Calculo de Kc
   ti=(k*kp)/(P2*tao);     //Calculo de ti
   td=0;
   
   //*************************************************************************//
   //*****************   DISEÑO POR CANCELACIÓN DE POLOS    *******************//
   //*************************************************************************//
   /*
   TsMA=7.5;                  //Tiempo deseado en Lazo Cerrado 
   kp=(tao)/(TsMA*k);      //Calculo de Kc
   ti=tao;                  //Calculo de Ti (Igual a la constante de tiempo)
   td=0;
   */
   //*************************************************************************//
   //*****************   SINTONIA POR ZIEGLER y NICHOLS    *******************//
   //*************************************************************************//
   
//!   kp=(0.9*tao)/(k*theta); //Z-N Rapido
//!   kp=kp/5;                //Z-N Lento
//!   ti=3.33*theta;
//!   td=0;
   
   //*************************************************************************//
   
  // Calculo do controle PID digital
   q0=kp*(1+T/(2*ti)+td/T);
   q1=-kp*(1-T/(2*ti)+(2*td)/T);
   q2=(kp*td)/T;
   
   //Inicializa el vector de control
   for(i=0;i<=kT;i++){
      u[i]=0;
   }
   
   //Inicializa las variables que contienen la lectura de salida
   yM=do_everything();
   yE=yM;yE_1=yM;yR=yM;yR_1=yM;
   while(1){ 
      yM=do_everything();
      //sprintf(msg,"%3.2f%cC\r\n",do_everything(),0xB0); 
      sprintf(msg,"I%3.2fFI%3.2fFC%ldRC%ldR",yM,yM,control/10,control/10); 
      printf("%s",msg); 
     // printf(" Hola Mundo ");
     
     if(muest>6){
       muest=0;
 
      //Salida del modelo nominal Pn(z)
      yE= b0*u[kT-d]+b1*u[kT-1-d]-a1*yE_1;
      
      //Diferencia entre proceso real y proceso nominal (error de prediccion)
      ep=yM-yE;
      
      //Salida del modelo nominal Rapido (sin retardo)
      yR= b0*u[kT]+b1*u[kT-1]-a1*yR_1;
      
      //Suma entre error de prediccion y salida predicha
      yp=yR+ep;
      
      //Llama la funcion del controlador PID
      PID();
      
      
      //desplaza el vector de la ley de control
      for(i=1;i<=kT;i++){
         u[i-1]=u[i];
      }
      
      //Actualizo los valores pasados
      yE_1=yE;
      yR_1=yR;
     }
      
      //tiempo de muestreo
      delay_ms(1000);
      muest++;
         
   } 
}

Código en Matlab

%% Ejemplo Monitoreo de señales en tiempo Real
function varargout=monitoreo(varargin)
parar=false;
fclose('all')
global tiempo salida escalon control
fig(1)=figure('name','Monitor','menubar','none','position',[200 200 800 700],'color',[0.9 0.6 0.3])
movegui(fig(1),'center');
axe(1)=axes('parent',fig(1),'units','pixels','position',[60 380 600 280],'xlim',[0 40],'ylim',[0 200],'xgrid','on','ygrid','on')
axe(2)=axes('parent',fig(1),'units','pixels','position',[60 50 600 280],'xlim',[0 40],'ylim',[0 100],'xgrid','on','ygrid','on')

set(get(axe(1),'XLabel'),'String','Tiempo (Seg)')
set(get(axe(1),'YLabel'),'String','Temperatura (°C)')

set(get(axe(2),'XLabel'),'String','Tiempo (Seg)')
set(get(axe(2),'YLabel'),'String','Control (%)')

lin(1)=line('parent',axe(1),'xdata',[],'ydata',[],'Color','r','LineWidth',2.5);
lin(2)=line('parent',axe(1),'xdata',[],'ydata',[],'Color','k','LineWidth',2);
lin(3)=line('parent',axe(2),'xdata',[],'ydata',[],'Color','r','LineWidth',2.5);


bot(1)=uicontrol('parent',fig(1),'style','pushbutton','string','Detener','position',[680 50 100 50],'callback',@stop,'fontsize',11)
bot(2)=uicontrol('parent',fig(1),'style','pushbutton','string','Enviar','position',[680 200 100 50],'callback',@enviar,'fontsize',11)

txbx(1)=uicontrol('parent',fig(1),'style','tex','string','Temp','position',[680 100 100 50],'fontsize',11)
txbx(2)=uicontrol('parent',fig(1),'style','edit','string','000','position',[680 250 100 50],'fontsize',11)



%% Funcion Pare
    function varargout=stop(hObject,evendata)
        parar=true;
        fclose(SerialP);
        delete(SerialP);
        clear SerialP;
        
    end

%% Funcion enviar
    function varargout=enviar(hObject,evendata)
        deg1=get(txbx(2),'string');
        deg=["S"+deg1+"$"];
        fwrite(SerialP,deg,'uchar');
    end
    %% funcion Graficar
   % function varargout=grafique(hObject,evendata)
        tiempo=[0];
        salida=[0];
        escalon=[0];
        control=[0];
        deg1="0";
   
        dt=1;
        limx=[0 40];
        limy=[0 200];
        set(axe(1),'xlim',limx,'ylim',limy);
        
        
    %% Configura el Puerto Serial
    
    SerialP=serial('COM8');
    set(SerialP,'Baudrate',9600); % se configura la velocidad a 9600 Baudios
    set(SerialP,'StopBits',1); % se configura bit de parada a uno
    set(SerialP,'DataBits',8); % se configura que el dato es de 8 bits, debe estar entre 5 y 8
    set(SerialP,'Parity','none'); % se configura sin paridad

    fopen(SerialP);
              
     %% Grafico
     k=5;nit = 10000;
        while(~parar)
%             if get(bot(3),'value')
%                drawnow(); % necesario para que actualice el grafico
%                continue % forza a salir del while para  la proxima interacion y de esta forma nop actualizo datos
%             
%             end


            % Lectura del Dato por Puerto Serial
            variable= (fread(SerialP,30,'uchar'));
            ini=find(variable==73); %Busca el I (Primer dato)
            ini=ini(1)+1;
            fin=find(variable==70); %Busca F (ultimo dato)
            fin= fin(find(fin>ini))-1;
            fin=fin(1);
            tempC=char(variable(ini:fin))';
            temp=str2num(tempC);
            
            %Lectura de la senal de control
            ini=find(variable==67); %Busca el C (Primer dato)
            ini=ini(1)+1;
            fin=find(variable==82); %Busca R (ultimo dato)
            fin= fin(find(fin>ini))-1;
            fin=fin(1);
            Con1=char(variable(ini:fin))';
            cont=str2num(Con1);
            
            
            set(txbx(1),'string',tempC);
            %Actualiza las variables del grafico
           
             tiempo=[tiempo tiempo(end)+dt];
             salida=[salida temp];
             control=[control cont];
             escalon=[escalon str2num(deg1)];
             set(lin(1),'xdata',tiempo,'ydata',salida);
             set(lin(2),'xdata',tiempo,'ydata',escalon);
             set(lin(3),'xdata',tiempo,'ydata',control);
             pause(dt); %% espera 1 seg para cada interación
             if tiempo(end)>=limx % actualizo grafica cuando llega a su limite en tiempo real
             limx=[0 limx(2)+40];
             set(axe(1),'xlim',limx) ;
             set(axe(2),'xlim',limx); 
             end
             
             if salida(end)>=limy % actualizo grafica cuando llega a su limite en tiempo real
             limy=[0 limy(2)+30];
             set(axe(1),'ylim',limy); 
             end
             
             k=k+1;
             if(k==nit)
                 parar=true;
             end
        end
        parar=false;
        

   % end
end

Perfecto controleros y controleras con esto llegamos al final de otro post del sitio web.

Recuerda que si deseas apoyar mi trabajo, puedes invitarme a un café y seguirme ayudando a mantener los servidores de este sitio web, es muy barato el café y contribuyes con el tiempo y esfuerzo invertidos en las clases elaboradas en el canal y pagina web: 👉Invitar a un Café a Sergio ☕️

Espero que esten muy bien y nos vemos en la próxima.

Volver al Curso de Microcontroladores PIC

Entradas relacionadas

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.

Comentarios (10)

Hola Sergio,

Agradezco que siempre estés compartiendo tu conocimiento. Tengo una duda, estoy realizando un sistema de control, atreves de un plc, pero la logica la estoy implementando por Simulink. Mi duda es: Ambas estructuras se encuentran bien hechas? la primera es para una simulacion, y la segunda, es la logica de control que estoy usando para mi proceso en tiempo real. Donde res calefactora, es el actuador y TT el sensor. Subnivel de temperatura es el error. y PV2(C) es la variable de proceso.

https://drive.google.com/file/d/1bcZd2G1xub_GLTUHAGsciBL5ELgIivc8/view?usp=sharing

Responder

Hola Christian, pues aparentemente parece que esta bien. Solo en la primera estructura modificas un poco la forma de presentarlo, el bloque predictor parece ser el mismo, ahora parece que usas un controlador en Serie, a diferencia del de abajo que es un controlador pid convensional.

Responder

Buenas noches ingeniero, primero felicitarle y agradecer por los conocimientos que nos proporciona, son de mucha ayuda para nosotros los que nos apasiona las leyes de control.
Tengo unas dudas, primero por que no usa un detector de cruce por cero para detectar el angulo de disparo y asi poder controlar (vi en otras referncias que usan ese tipo de metodo para controlar cargas AC). Existe alguna diferencia sustancial que en este caso solo usemos rele de estado solido que recibe la señal de PWM del pic para controlar la temperatura. No es necesario el cruce por cero???
Segundo si yo tengo un crystal de 48MHz la minima frecuncia que tengo es 2031Hz para PWM y como podemos observar es mucha mas grande a comparacion del la fercuencia q nos suministra la red electrica, haciendo unas pruebas en un sistema de que es una resistencia (tipo plancha domestica) y un potenciometro para cambiar la temperatura no estoy consiguiendo un cambio.
Estoy haciendo algo mal?
De antemano muchas gracias por la respuesta, saludos desde La Paz Bolivia.

Responder

Hola Juan, mi recomendación es que uses el detector por cruce por cero, tendrás un mejor control, eso es un proyecto que voy a hacer en un futuro en el canal, de hecho ya tengo todo comprado para hacerlo. El dispositivo que usé en el relé de estado solido (el cual compré hace muchos años) internamente tiene implementado un detector de cruce por cero según el fabricante y como fue un proyecto que había hecho hace tiempo y funcionó entonces lo traje al canal. Por eso es que no lo detallé aquí, mas adelante pretendo hablar más sobre eso.
No se como tengas configurado tu PWM, tienes que fijarte en la resolución que tienes según tu configuración para asignar un valor adecuado al duty cicle. Por ejemplo el PWM de esta entrada fue de 1000hz (con cristal de 4MHz), por lo tanto la resolución va de 0 – 1000, que son los valores que puedo colocar en el duty cicle.
En la entrada de PWM, explica eso y también muestra un ejemplo donde se configura el pwm a 4000hz donde la resolución se reduce hasta 63. Es decir que en el duty cicle solo puedo enviar valores de 0 a 63. Saludos.

Responder

Boa tarde, o que faz a função do_everything? Não tem o protótipo da função.

Responder

Olá Glaucio ela lê a temperatura da termocupla, mas isso está explicado na entrada referente à leitura do sensor no seguinte link.

Responder

Un saludo Sergio Castaño. Gracias por compartir tus conocimientos. Y así capcitarnos en al área de control y programación.

Tengo una duda con respecto a este ciclo for…

for(i=0;i<5;i++){
if(dat[i]=='S'){
ini=i+1;
i=10;
}
}

Mi pregunta es… ¿ por que al final pones que i=10; ?.

Responder

Hola Manuel, al final coloco i=10 es para que se salga del for. pues en el arreglo dat estoy buscando S o $, cuando los encuentro, no me interesa seguir buscando, por eso coloco i=10 para salir del for. Saludos.

Responder

tus videos me encantan, eh aprendido mucho con ellos, pero en estos últimos videos si que me he perdido, quisiera implentar el control pid solo con la termocupla y la resistencia sin la parte de matlab pero no me doy idea de como si pudieras ayudarme te estaría muy agradecido, de antemano gracias

Responder

Juan es exactamente igual, en tu caso solo vas a usar la resistencia, el relé de estado solido y la termocupla. Omites la parte de MATLAB, que únicamente es para visualización. En tu caso puedes usar un LCD para ver los resultados, de la forma como lo hicimos en los videos de los PID anteriores.

Responder