Pdepe solver be run in loop for designing a control scheme.

Hi, i am designing a control scheme for a thermal system defined by a 2nd order pde. I used pdepe solver to get the solution but driven by input voltage. But in-order-to design a control scehme i need to run the pdepe solver in iterations .please guide.

[Torsten]

What are the difficulties you encounter or expect to encounter ?

[maria atiq]

When i try to give index to the sol() term . Matlab gives the error that it is illegal to index this function.

[Torsten]

You can save "sol" as a cell array in each iteration as "sol{i}".

But why do you want to save intermediate sol terms if you iterate and get a new "sol" in the next iteration ? Do you want to have a kind of convergence monitor after exiting the iteration loop ?

[maria atiq]

No, actually my equation has voltage as an input control term and the PID i am implementing via code , Updates that voltage according to the error in between T out put from sol() and T desired.

But the PID code needs to run the sol() multiple times with every iteration of voltage , until the error becomes zero

[Torsten]

I think we are unable to help you without some code that clarifies your problem.

[maria atiq]

I am using a simple PID code , here it is used for a dummy heat equation T(i+1). where as i want to run it for my 2nd order pde which i am solving seperately using a pdepe solver.

for i = 1:n
    Error(i+1) = T_desired - T(i); % error entering the PID controller
    
    Prop(i+1) = Error(i+1);% error of proportional term
    Der(i+1)  = (Error(i+1) - Error(i))/dt; % derivative of the error
    Int(i+1)  = (Error(i+1) + Error(i))*dt/2; % integration of the error
    I(i+1)    = sum(Int); % the sum of the integration of the error
    
    PID(i+1)  = Kp*Prop(i) + Ki*I(i+1)+ Kd*Der(i); % the three PID terms
    
    
    STATE1(i+1) = sum(PID); % sum PID term to calculate the first integration
    voltage(i+1) = (STATE1(i+1) + STATE1(i))*dt/2; 

    T(i+1) = ((voltage(i+1)^2)/(R))*(1-alpha*(T(i)-Tref));
    
end

[maria atiq]

I was thinking to write it as

sol(i+1) = pdepe(equfn, icfn, bcfn...) T(i+1)=sol(:,:,1);

Eqfun takes in the iterative voltage values

[Torsten]

What does the index i stand for ? The average solution for T over all x-values at time (i*dt) ?

[maria atiq]

the solution of T at said input voltage(i)

[maria atiq]

i am getting the followimg error when i try to index sol

Unable to perform assignment because the left and right sides have a different number of elements.

Error in transistor_example (line 42)

sol(i+1) = pdepe(m,eqn,ic,@transistorBC,x,t);

[maria atiq]

it is in the eqnfcn called in pdepe

sol=pdepe(eqnfcn---)

[Torsten]

Again my question: What does the index i stand for ? If you solve the heat equation with "pdepe", you get a matrix T as solution where the rows of T stand for the time points t and the columns stand for grid points x. You only use a single scalar T(i+1). I cannot make sense of it.

[maria atiq]

Yes i agree that sol gives T values against t and x.

I want to run the solver in loop for varying input voltages that are going to be input in eqn function of solver . So for a series of voltages or let say updating voltages the solver should run on its own

[Torsten]

You still haven't told us what the i stands for. Is it the ith x-value, the ith time value, the ith iteration ? And how does your computed voltage(i) influence the computed T from "pdepe" ?

[maria atiq]

Hi, i am attaching my code.

dt = 0.000162;          % sampling time
Time = 0.01;          % total simulation time in seconds
n = round(Time/dt);
prompt="What is the desired temp?"
T_desired = input(prompt);        % desired output, or reference point

Kp = 1;             % proportional term Kp
Ki = 0.5005;           % Integral term Ki
Kd = 0.001;           % derivative term Kd


% pre-assign all the arrays to optimize simulation time
Prop(1:k+1) = 0; Der(1:k+1) = 0; Int(1:k+1) = 0; I(1:k+1) = 0;
PID(1:k+1) = 0;
Temp_avg(1:k+1) = 0;
voltage (1:k+1) = 0;
Error(1:k+1) = 0;
state1(1:k+1) = 0; STATE1(1:k+1) = 0;
state2(1:k+1) = 0; STATE2(1:k+1) = 0;

for i=1:k
 

   Error(i)=T_desired- Temp_avg(i);
Prop(i+1) = Error(i+1);% error of proportional term
    Der(i+1)  = (Error(i+1) - Error(i))/dt; % derivative of the error
    Int(i+1)  = (Error(i+1) + Error(i))*dt/2; % integration of the error
    I(i+1)    = sum(Int); % the sum of the integration of the error
    
    PID(i+1)  = Kp*Prop(i) + Ki*I(i+1)+ Kd*Der(i); % the three PID terms
    
    
    STATE1(i+1) = sum(PID); % sum PID term to calculate the first integration
    state2(i+1) = (STATE1(i+1) + STATE1(i))*dt/2; % output after the first integrator
    STATE2(i+1) = sum(state2); % sum output of first integrator to calculate the second integration
    voltage(i+1) = (STATE2(i+1) + STATE2(i))*dt/2; 
C.V=voltage(i+1);
x1=linspace(0,55,10);
x2=linspace(55.01,59,100);
x3=linspace(59.01,589,100);
x=cat(2,x1,x2,x3); 
t = linspace(0,0.01,60);
m = 1;
eqn = @(x,t,T,dTdx) heateqfull1(x,t,T,dTdx,C);
ic = @(x) incondfull1(x);
sol = pdepe(m,eqn,ic,@boundarycondfull1,x,t);
Tk=sol(:,:,1);
temp=Tk-273.15; 
Temp_heater=temp(1:60,10:110);
Temp_avg=mean(Temp_heater,2);

end

[maria atiq]

I want my pdepe solver to run in loop for every calculated value of voltage untill error is reduced to zero. but MATLAB does not allow me to index sol function

[Torsten]

The code is incomplete and cannot be executed.

The variable k is not defined and a reasonable value for T_desired is needed. Further,the functions for pdepe are missing.

I suggest testing the code before posting it.

[maria atiq]

dt = 0.000168;          % sampling time
Time = 0.01;          % total simulation time in seconds
k = round(Time/dt);

T_desired = 200;        % desired output, or reference point

Kp = 1;             % proportional term Kp
Ki = 0.5005;           % Integral term Ki
Kd = 0.001;           % derivative term Kd


% pre-assign all the arrays to optimize simulation time
Prop(1:k+1) = 0; Der(1:k+1) = 0; Int(1:k+1) = 0; I(1:k+1) = 0;
PID(1:k+1) = 0;
Temp_avg(1:k+1) = 0;
voltage (1:k+1) = 0;
Error(1:k+1) = 0;
state1(1:k+1) = 0; STATE1(1:k+1) = 0;
state2(1:k+1) = 0; STATE2(1:k+1) = 0;

for i=1:k
 
%    T_desired=145;
   Error(i)=T_desired- Temp_avg(i);
%     error(i)=abs(Error(i)/100);
Prop(i+1) = Error(i+1);% error of proportional term
    Der(i+1)  = (Error(i+1) - Error(i))/dt; % derivative of the error
    Int(i+1)  = (Error(i+1) + Error(i))*dt/2; % integration of the error
    I(i+1)    = sum(Int); % the sum of the integration of the error
    
    PID(i+1)  = Kp*Prop(i) + Ki*I(i+1)+ Kd*Der(i); % the three PID terms
    
    %% You can replace the follwoing five lines with your system/hardware/model
    STATE1(i+1) = sum(PID); % sum PID term to calculate the first integration
    state2(i+1) = (STATE1(i+1) + STATE1(i))*dt/2; % output after the first integrator
    STATE2(i+1) = sum(state2); % sum output of first integrator to calculate the second integration
    voltage(i+1) = (STATE2(i+1) + STATE2(i))*dt/2; 
C.V=voltage(i+1);
x1=linspace(0,55,10);
x2=linspace(55.01,59,100);
x3=linspace(59.01,589,100);
x=cat(2,x1,x2,x3); 
t = linspace(0,0.01,60);
m = 1;
eqn = @(x,t,T,dTdx) heateqfull1(x,t,T,dTdx,C);
ic = @(x) incondfull1(x);
sol = pdepe(m,eqn,ic,@boundarycondfull1,x,t);
Tk=sol(:,:,1);
temp=Tk-273.15; 
Temp_heater=temp(1:60,10:110);
Temp_avg=mean(Temp_heater,2);

end
T = 0:dt:Time;
Reference = T_desired*ones(1,i);
plot(T,Reference,'r',T,Temp_avg,'b');
xlabel('Time (sec)')
legend('Desired','Simulated')

function [c,f,s] = heateqfull1(x,t,T,dTdx,C)
           %%---------- universal parametrs -------------
tm= 1.7;                            % thickness of membrane in micrometer 
th= 0.3;                            % thickness of heater in micrometer
rh= 55;                             % radius of heater in micrometer
rm= 589;                            % radius of membrane in micrometer 
wh= 4;                              % width of heater in micrometer 
r1= 55;                             % radius of region 1 in micrometer 
r2= 59;                             % radius of region 2 in micrometer 
r3= 589;                            % radius of region 3 in micrometers
k= 1.61*10^-5;                      % effective thermal conductivity of Si3N4 and SiO2 in W/µmK 
ka= 2.623 * 10^-8;                  % thermal conductivity of air in W/µmK 
kw= 0.598 * 10^-6;                  % thermal conductivity of air in W/µmK 
hc= 2*10^-10;                       % effective heat trasnfer co-efficent in W/µm2K
rho= 2.76*10^-12;                   % denisty of membrane g/um3
shc= 0.68812;                       % specific heat capacity of membrane in J/(g*K)
Ta= 293.15;                         % room temperature in K 
kpt= 2.95*10^-5;                    % effective thermal conductivity of platinum in W/µmK

     %%__ resistance of heater depending upon thickness and resistivity ___
         
rrho= 0.105;                        % resistivity of platinum in ohm um ;  1.05*10*-7 ohm m 2014 paper 
St=1;                               %2013
Sc=1;                               %2013
Wc=24;                              % extrernal wide width contact (nc-ext *wti=12*2=24 )2013
Sj=20;                              % length of external wide width contact (ns-ext *wti=10*2=20)2013
Sa=6;                               % 2015 paper supporting documnet as this parameter is not in 2013 paper 
Lc=503;                             % 2015 paper supporting documnet as this parameter is not in 2013 paper
th=0.3;                             % thickness of heater in um from 2015 paper 
Ro=rrho*([(pi*((5*rh)+(6*St)+((25/4)*wh)-(4*Sc)+(2*Sj))/(wh*th))]+[(2*Lc)/(Wc*th)]);
alpha = 3.927*10^-3;                % temperature coefficient of resistivity in 1/K
V =C.V;                           % voltage for input in Volts 

         %%--------  specific user defined parameters-----------
ci= (rho*shc)/k;            
si=(2*hc)/(k*tm);
sii=1/(2*kpt*pi*r1*wh*tm*Ro);

c = ci;
if x>=0 & x<=55
         s=-si *(T-293.15);
         f = dTdx;
      else if x>=55 & x<=59
         s = -(si*(T-293.15))+ [V^2 *(1-(alpha*(T-293.15)))]/sii ;
         f = dTdx;
      else if x>=59
              s=-si *(T-293.15);
              f = dTdx;
          end
          end
end 
end
function [pl,ql,pr,qr] = boundarycondfull1(xl,Tl,xr,Tr,t)

  pl = 0;
  ql = 1;
  pr = Tr-293.15  ;
  qr = 0;


end
function T0 = incondfull1(x)

T0 = 293.15;

end

[maria atiq]

Hey, can you run it now?

[Torsten]

Ok, I see that your PDE depends somehow on the input C, and you say you want to change C in the loop. But how ?

[maria atiq]

Basicaly, C is a structure aray that is used enter value into the solver. "create a structure array with fields for each phyical parameter. When you later define functions for the equations, the initial condition, and the boundary conditions, you can pass in this structure as an extra argument so that the functions have access to the parameter."

So, the updated voltage value is saved as C.V and it is callled in equation function.

[maria atiq]

the function is called in the solver .

[Torsten]

So, the updated voltage value is saved as C.V and it is callled in equation function.

But you don't update C.V depending on sol, do you ?

[maria atiq]

ideally the C.V should be updated as the alue of error reduces to zero and for every iteration we have a new value of C.V .

here is the problem, this value of C.V should be used in new sol at every ieration to generate an updated Temp_avg. untill error is zero. and desired temp is achieved.

[Torsten]

I think this is somehow what you want.

I tried to keep "voltage" of size k, but the PID control as I interpreted it from your code gives me an array of size k-2 in the end which leads to a size dimension mismatch in function "heateqfull1".

Hope you can solve the problem.

dt = 0.000168;          % sampling time
Time = 0.01;          % total simulation time in seconds
k = round(Time/dt);
T_desired = 200;        % desired output, or reference point
Kp = 1;             % proportional term Kp
Ki = 0.5005;           % Integral term Ki
Kd = 0.001;           % derivative term Kd
% pre-assign all the arrays to optimize simulation time
Prop(1:k+1) = 0; Der(1:k+1) = 0; Int(1:k+1) = 0; I(1:k+1) = 0;
PID(1:k+1) = 0;
Temp_avg(1:k+1) = 0;
voltage (1:k) = 0;
Error(1:k+1) = Inf;
state1(1:k+1) = 0; STATE1(1:k+1) = 0;
state2(1:k+1) = 0; STATE2(1:k+1) = 0;
Times = linspace(0,Time,k+1);
m = 1;
x1 = linspace(0,55,10);
x2 = linspace(55.01,59,100);
x3 = linspace(59.01,589,100);
x = cat(2,x1,x2,x3);
ic = @(x) incondfull1(x);
size(Times)

ans = 1×2

1 61

size(voltage)

ans = 1×2

1 60

while max(abs(Error))>1e-2
   eqn = @(x,t,T,dTdx) heateqfull1(x,t,T,dTdx,Times,voltage);
   m = 1;   
   sol = pdepe(m,eqn,ic,@boundarycondfull1,x,Times);
   for i = 1:k+1
     Temp_avg(i) = 2/589^2*trapz(sol(i,:,1).*x);
   end
   Error = T_desired - Temp_avg;
   Prop = Error(2:end);
   Der = (Error(2:end)-Error(1:end-1))/dt;
   Int = (Error(2:end)-Error(1:end-1))*dt/2;
   I = cumsum(Int);
   PID = Kp*Prop + Ki*I + Kd*Der;
   STATE1 = cumsum(PID);
   state2 = (STATE1(2:end)+STATE1(1:end-1))*dt/2;
   STATE2 = cumsum(state2);
   voltage = (STATE2(2:end) + STATE2(1:end-1))*dt/2;
   size(voltage)
end

ans = 1×2

1 58

Index exceeds the number of array elements. Index must not exceed 58.

Error in solution>heateqfull1 (line 88)
V = voltage_inter(j);

Error in solution>@(x,t,T,dTdx)heateqfull1(x,t,T,dTdx,Times,voltage) (line 26)
eqn = @(x,t,T,dTdx) heateqfull1(x,t,T,dTdx,Times,voltage);

Error in pdepe/pdeodes (line 353)
[cR,fR,sR] = feval(pde,xi(ii),tnow,U,Ux,varargin{:});

Error in ode15s (line 553)
rhs = hinvGak*ode(tnew,ynew) - Mtnew*(psi+difkp1);

Error in pdepe (line 287)
[t,y] = ode15s(@pdeodes,t,y0(:),opts);

function [c,f,s] = heateqfull1(x,t,T,dTdx,t_inter,voltage_inter)
           %%---------- universal parametrs -------------
tm= 1.7;                            % thickness of membrane in micrometer 
th= 0.3;                            % thickness of heater in micrometer
rh= 55;                             % radius of heater in micrometer
rm= 589;                            % radius of membrane in micrometer 
wh= 4;                              % width of heater in micrometer 
r1= 55;                             % radius of region 1 in micrometer 
r2= 59;                             % radius of region 2 in micrometer 
r3= 589;                            % radius of region 3 in micrometers
k= 1.61*10^-5;                      % effective thermal conductivity of Si3N4 and SiO2 in W/µmK 
ka= 2.623 * 10^-8;                  % thermal conductivity of air in W/µmK 
kw= 0.598 * 10^-6;                  % thermal conductivity of air in W/µmK 
hc= 2*10^-10;                       % effective heat trasnfer co-efficent in W/µm2K
rho= 2.76*10^-12;                   % denisty of membrane g/um3
shc= 0.68812;                       % specific heat capacity of membrane in J/(g*K)
Ta= 293.15;                         % room temperature in K 
kpt= 2.95*10^-5;                    % effective thermal conductivity of platinum in W/µmK
     %%__ resistance of heater depending upon thickness and resistivity ___
         
rrho= 0.105;                        % resistivity of platinum in ohm um ;  1.05*10*-7 ohm m 2014 paper 
St=1;                               %2013
Sc=1;                               %2013
Wc=24;                              % extrernal wide width contact (nc-ext *wti=12*2=24 )2013
Sj=20;                              % length of external wide width contact (ns-ext *wti=10*2=20)2013
Sa=6;                               % 2015 paper supporting documnet as this parameter is not in 2013 paper 
Lc=503;                             % 2015 paper supporting documnet as this parameter is not in 2013 paper
th=0.3;                             % thickness of heater in um from 2015 paper 
Ro=rrho*([(pi*((5*rh)+(6*St)+((25/4)*wh)-(4*Sc)+(2*Sj))/(wh*th))]+[(2*Lc)/(Wc*th)]);
alpha = 3.927*10^-3;                % temperature coefficient of resistivity in 1/K
%V =C.V;                           % voltage for input in Volts 
         %%--------  specific user defined parameters-----------
ci= (rho*shc)/k;            
si=(2*hc)/(k*tm);
sii=1/(2*kpt*pi*r1*wh*tm*Ro);
c = ci;
if x>=0 & x<=55
         s=-si *(T-293.15);
         f = dTdx;
      else if x>=55 & x<=59
         for j = 1:numel(t_inter)-1
           if t >= t_inter(j) & t <= t_inter(j+1)
             V = voltage_inter(j);
             break
           end
         end
         s = -(si*(T-293.15))+ [V^2 *(1-(alpha*(T-293.15)))]/sii ;
         f = dTdx;
      else if x>=59
              s=-si *(T-293.15);
              f = dTdx;
          end
          end
end 
end
function [pl,ql,pr,qr] = boundarycondfull1(xl,Tl,xr,Tr,t)
  pl = 0;
  ql = 1;
  pr = Tr-293.15  ;
  qr = 0;
end
function T0 = incondfull1(x)
T0 = 293.15;
end

[maria atiq]

Hey, thank you for your help. I sorted the size issue, however still getting the index error.

where in heateqfull1 V=voltage_inter(j)

it gives the following error "unable to perform proper assignment"

clc
clear all

dt = 0.000168;          % sampling time
Time = 0.01;          % total simulation time in seconds
k = round(Time/dt);
T_desired = 200;        % desired output, or reference point
Kp = 1;             % proportional term Kp
Ki = 0.5005;           % Integral term Ki
Kd = 0.001;           % derivative term Kd
% pre-assign all the arrays to optimize simulation time
Prop(1:k+1) = 0; Der(1:k+1) = 0; Int(1:k+1) = 0; I(1:k+1) = 0;
PID(1:k+1) = 0;
Temp_avg(1:k+1) = 0;
voltage (1:k) = 0;
Error(1:k+1) = Inf;
state1(1:k+1) = 0; STATE1(1:k+1) = 0;
state2(1:k+1) = 0; STATE2(1:k+1) = 0;
% Times = linspace(0,Time,k+1);
t = linspace(0,0.01,60);
m = 1;
x1 = linspace(0,55,10);
x2 = linspace(55.01,59,100);
x3 = linspace(59.01,589,100);
x = cat(2,x1,x2,x3);
ic = @(x) incondfull1(x);
size(t)

ans = 1×2

1 60

% ans = 1×2
%      1    61
size(voltage)

ans = 1×2

1 60

% ans = 1×2
%      1    60
while max(abs(Error))>1e-2
   eqn = @(x,t,T,dTdx) heateqfull1(x,t,T,dTdx,t,voltage);
   m = 1;   
   sol = pdepe(m,eqn,ic,@boundarycondfull1,x,t);
   for i = 1:k+1
     Temp_avg(i) = 2/589^2*trapz(sol(i,:,1).*x);
   end
   Error = T_desired - Temp_avg;
   Prop = Error(2:end);
   Der = (Error(2:end)-Error(1:end-1))/dt;
   Int = (Error(2:end)-Error(1:end-1))*dt/2;
   I = cumsum(Int);
   PID = Kp*Prop + Ki*I + Kd*Der;
   STATE1 = cumsum(PID);
   state2 = (STATE1(2:end)+STATE1(1:end-1))*dt/2;
   STATE2 = cumsum(state2);
   voltage = (STATE2(2:end) + STATE2(1:end-1))*dt/2;
   size(voltage)
end

Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 0-by-0.

Error in pdepe/pdeodes (line 357)
up(:,ii) = ((xim(ii) * fR - xim(ii-1) * fL) + ...

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 148)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in pdepe (line 287)
[t,y] = ode15s(@pdeodes,t,y0(:),opts);

% ans = 1×2
%      1    58
% Index exceeds the number of array elements. Index must not exceed 58.
% 
% Error in solution>heateqfull1 (line 88)
%              V = voltage_inter(j);
% 
% Error in solution>@(x,t,T,dTdx)heateqfull1(x,t,T,dTdx,Times,voltage) (line 26)
%    eqn = @(x,t,T,dTdx) heateqfull1(x,t,T,dTdx,Times,voltage);
% 
% Error in pdepe/pdeodes (line 353)
%             [cR,fR,sR] = feval(pde,xi(ii),tnow,U,Ux,varargin{:});
% 
% Error in ode15s (line 553)
%                 rhs = hinvGak*ode(tnew,ynew) -  Mtnew*(psi+difkp1);

% Error in pdepe (line 287)
%         [t,y] = ode15s(@pdeodes,t,y0(:),opts);
function [c,f,s] = heateqfull1(x,t,T,dTdx,t_inter,voltage_inter)
           %%---------- universal parametrs -------------
tm= 1.7;                            % thickness of membrane in micrometer 
th= 0.3;                            % thickness of heater in micrometer
rh= 55;                             % radius of heater in micrometer
rm= 589;                            % radius of membrane in micrometer 
wh= 4;                              % width of heater in micrometer 
r1= 55;                             % radius of region 1 in micrometer 
r2= 59;                             % radius of region 2 in micrometer 
r3= 589;                            % radius of region 3 in micrometers
k= 1.61*10^-5;                      % effective thermal conductivity of Si3N4 and SiO2 in W/µmK 
ka= 2.623 * 10^-8;                  % thermal conductivity of air in W/µmK 
kw= 0.598 * 10^-6;                  % thermal conductivity of air in W/µmK 
hc= 2*10^-10;                       % effective heat trasnfer co-efficent in W/µm2K
rho= 2.76*10^-12;                   % denisty of membrane g/um3
shc= 0.68812;                       % specific heat capacity of membrane in J/(g*K)
Ta= 293.15;                         % room temperature in K 
kpt= 2.95*10^-5;                    % effective thermal conductivity of platinum in W/µmK
     %%__ resistance of heater depending upon thickness and resistivity ___
         
rrho= 0.105;                        % resistivity of platinum in ohm um ;  1.05*10*-7 ohm m 2014 paper 
St=1;                               %2013
Sc=1;                               %2013
Wc=24;                              % extrernal wide width contact (nc-ext *wti=12*2=24 )2013
Sj=20;                              % length of external wide width contact (ns-ext *wti=10*2=20)2013
Sa=6;                               % 2015 paper supporting documnet as this parameter is not in 2013 paper 
Lc=503;                             % 2015 paper supporting documnet as this parameter is not in 2013 paper
th=0.3;                             % thickness of heater in um from 2015 paper 
Ro=rrho*([(pi*((5*rh)+(6*St)+((25/4)*wh)-(4*Sc)+(2*Sj))/(wh*th))]+[(2*Lc)/(Wc*th)]);
alpha = 3.927*10^-3;                % temperature coefficient of resistivity in 1/K
%V =C.V;                           % voltage for input in Volts 
         %%--------  specific user defined parameters-----------
ci= (rho*shc)/k;            
si=(2*hc)/(k*tm);
sii=1/(2*kpt*pi*r1*wh*tm*Ro);
c = ci;
if x>=0 & x<=55
         s=-si *(T-293.15);
         f = dTdx;
      else if x>=55 & x<=59
         for j = 1:numel(t_inter)-1
           if t >= t_inter(j) & t <= t_inter(j+1)
             V = voltage_inter(j);
             break
           end
         end
         s = -(si*(T-293.15))+ [voltage_inter(j).^2 *(1-(alpha*(T-293.15)))]/sii ;
         f = dTdx;
      else if x>=59
              s=-si *(T-293.15);
              f = dTdx;
          end
          end
end 
end
function [pl,ql,pr,qr] = boundarycondfull1(xl,Tl,xr,Tr,t)
  pl = 0;
  ql = 1;
  pr = Tr-293.15  ;
  qr = 0;
end
function T0 = incondfull1(x)
T0 = 293.15;
end

[Torsten]

I can't understand why you changed

Times = linspace(0,Time,k+1);

to

t = linspace(0,0.01,60);

It makes no sense at all for me.

Given time points t0,...,tk and average temperatures T0,...,Tk over the domain at these time points, you need to compute voltages V1,...,Vk for the time intervals (V1 applied in [t0 t1], V2 applied in [t1 t2],...,Vk applied in [t(k-1) tk]) by your optimal control equations. I made an attempt to get k voltages from your equations, but only got k-2. So you will need to check the commands in the while-loop.

[maria atiq]

Hey, agreed. I tried looking into it but couldnt sort the indexing issue.

SO, i tried adding extra zeros at the end of the voltage (size 58) just to check the rest of the code.

But it goes into the while loop and keeps running and doesnt come out of it. Can you check what is the issue?

dt = 0.000168;          % sampling time
Time = 0.01;          % total simulation time in seconds
k = round(Time/dt);
T_desired = 100;        % desired output, or reference point
Kp = 1;             % proportional term Kp
Ki = 0.5005;           % Integral term Ki
Kd = 0.001;           % derivative term Kd
% pre-assign all the arrays to optimize simulation time
Prop(1:k+1) = 0; Der(1:k+1) = 0; Int(1:k+1) = 0; I(1:k+1) = 0;
PID(1:k+1) = 0;
Temp_avg(1:k+1) = 0;
voltage (1:k) = 0;
VOLTAGE(1:k)=0;
Error(1:k+1) = Inf;
state1(1:k+1) = 0; STATE1(1:k+1) = 0;
state2(1:k+1) = 0; STATE2(1:k+1) = 0;
Times = linspace(0,Time,k+1);
m = 1;
x1 = linspace(0,55,10);
x2 = linspace(55.01,59,100);
x3 = linspace(59.01,589,100);
x = cat(2,x1,x2,x3);
ic = @(x) incondfull1(x);
size(Times)
size(voltage)
while max(abs(Error))>1e-2
   eqn = @(x,t,T,dTdx) heateqfull1(x,t,T,dTdx,Times,VOLTAGE);
   m = 1;   
   sol = pdepe(m,eqn,ic,@boundarycondfull1,x,Times);
   for i = 1:k+1
     Temp_avg(i) = 2/589^2*trapz(sol(i,:,1).*x);
   end
   Error = T_desired - Temp_avg;
   Prop = Error(2:end);
   Der = (Error(2:end)-Error(1:end-1))/dt;
   Int = (Error(2:end)-Error(1:end-1))*dt/2;
   I = cumsum(Int);
   PID = Kp*Prop + Ki*I + Kd*Der;
   STATE1 = cumsum(PID);
   state2 = (STATE1(2:end)+STATE1(1:end-1))*dt/2;
   STATE2 = cumsum(state2);
   voltage = (STATE2(2:end) + STATE2(1:end-1))*dt/2;
   size(voltage)
   VOLTAGE=[voltage 0 0]
   size(VOLTAGE)
end


function [c,f,s] = heateqfull1(x,t,T,dTdx,t_inter,voltage_inter)
           %%---------- universal parametrs -------------
tm= 1.7;                            % thickness of membrane in micrometer 
th= 0.3;                            % thickness of heater in micrometer
rh= 55;                             % radius of heater in micrometer
rm= 589;                            % radius of membrane in micrometer 
wh= 4;                              % width of heater in micrometer 
r1= 55;                             % radius of region 1 in micrometer 
r2= 59;                             % radius of region 2 in micrometer 
r3= 589;                            % radius of region 3 in micrometers
k= 1.61*10^-5;                      % effective thermal conductivity of Si3N4 and SiO2 in W/µmK 
ka= 2.623 * 10^-8;                  % thermal conductivity of air in W/µmK 
kw= 0.598 * 10^-6;                  % thermal conductivity of air in W/µmK 
hc= 2*10^-10;                       % effective heat trasnfer co-efficent in W/µm2K
rho= 2.76*10^-12;                   % denisty of membrane g/um3
shc= 0.68812;                       % specific heat capacity of membrane in J/(g*K)
Ta= 293.15;                         % room temperature in K 
kpt= 2.95*10^-5;                    % effective thermal conductivity of platinum in W/µmK
     %%__ resistance of heater depending upon thickness and resistivity ___
         
rrho= 0.105;                        % resistivity of platinum in ohm um ;  1.05*10*-7 ohm m 2014 paper 
St=1;                               %2013
Sc=1;                               %2013
Wc=24;                              % extrernal wide width contact (nc-ext *wti=12*2=24 )2013
Sj=20;                              % length of external wide width contact (ns-ext *wti=10*2=20)2013
Sa=6;                               % 2015 paper supporting documnet as this parameter is not in 2013 paper 
Lc=503;                             % 2015 paper supporting documnet as this parameter is not in 2013 paper
th=0.3;                             % thickness of heater in um from 2015 paper 
Ro=rrho*([(pi*((5*rh)+(6*St)+((25/4)*wh)-(4*Sc)+(2*Sj))/(wh*th))]+[(2*Lc)/(Wc*th)]);
alpha = 3.927*10^-3;                % temperature coefficient of resistivity in 1/K
%V =C.V;                           % voltage for input in Volts 
         %%--------  specific user defined parameters-----------
ci= (rho*shc)/k;            
si=(2*hc)/(k*tm);
sii=1/(2*kpt*pi*r1*wh*tm*Ro);
c = ci;
if x>=0 & x<=55
         s=-si *(T-293.15);
         f = dTdx;
      else if x>=55 & x<=59
         for j = 1:numel(t_inter)-1
           if t >= t_inter(j) & t <= t_inter(j+1)
             V = voltage_inter(j);
             break
           end
         end
         s = -(si*(T-293.15))+ [V^2 *(1-(alpha*(T-293.15)))]/sii ;
         f = dTdx;
      else if x>=59
              s=-si *(T-293.15);
              f = dTdx;
          end
          end
end 
end
function [pl,ql,pr,qr] = boundarycondfull1(xl,Tl,xr,Tr,t)
  pl = 0;
  ql = 1;
  pr = Tr-293.15  ;
  qr = 0;
end
function T0 = incondfull1(x)
T0 = 293.15;
end

[maria atiq]

It would mean alot if this issue can be sorted as i have been stuck at it for way too long. Kind regards.

[Torsten]

You must think about what you want to achieve in your voltage control: You want to find a time-dependent voltage applied for 55<=x<=59 such that nearly from the first time instant on, the average temperature over the length is constantly at T_desired = 100. For me it sounds that this is impossible to achieve.

[Bill Greene]

I know I'm late to this conversation but I was curious so ran the code.

When I print

max(abs(Error))

at the bottom of the while loop, it equals 100 for all iterations.

PS: Is there an introductory reference that describes the controller you are trying to implement?

[maria atiq]

I want to implement a very basic PID , but since systems is not known in the form of tf or SS. I can't use predefined commands for PID, hence i wrote the code for PID. I want to run the sol ( pdepe) solver in loop untill error is reduced to zero.

[Torsten]

I made an error in the code. It must read

Temp_avg(i) = 2/589^2*trapz(x,sol(i,:,1).*x);

instead of

Temp_avg(i) = 2/589^2*trapz(sol(i,:,1).*x);

I don't have an idea about control theory - so I don't know what's wrong about your PID control scheme.

[maria atiq]

The control part i have tested on a simple thermal equation and it works fine.

clc;clear;close all


T_desired = 180;        % desired output, or reference point


R=108.6830;
alpha = 3.927*10^-3; 
Tref=293.15;
V=0.5;
Kp = 1;             % proportional term Kp
Ki = 0.5;           % Integral term Ki
Kd = 0.01;           % derivative term Kd
dt = 0.01;          % sampling time
Time = 10;          % total simulation time in seconds
n = round(Time/dt); % number of samples

% pre-assign all the arrays to optimize simulation time
Prop(1:n+1) = 0; Der(1:n+1) = 0; Int(1:n+1) = 0; I(1:n+1) = 0;
PID(1:n+1) = 0;
Temp(1:n+1) = 0;
Output(1:n+1) = 0;
Error(1:n+1) = 0;
state1(1:n+1) = 0; STATE1(1:n+1) = 0;
state2(1:n+1) = 0; STATE2(1:n+1) = 0;

for i = 1:n
    Error(i+1) = T_desired - Temp(i); % error entering the PID controller
    
    Prop(i+1) = Error(i+1);% error of proportional term
    Der(i+1)  = (Error(i+1) - Error(i))/dt; % derivative of the error
    Int(i+1)  = (Error(i+1) + Error(i))*dt/2; % integration of the error
    I(i+1)    = sum(Int); % the sum of the integration of the error
    
    PID(i+1)  = Kp*Prop(i) + Ki*I(i+1)+ Kd*Der(i); % the three PID terms
    
    
    STATE1(i+1) = sum(PID); % sum PID term to calculate the first integration
    Output(i+1) = (STATE1(i+1) + STATE1(i))*dt/2; % output after the first integrator
%      
    Temp(i+1) = ((Output(i+1).^2)/R)*(1-alpha*(Temp(i)-Tref));
end



% plot results
T = 0:dt:Time;
Reference = T_desired*ones(1,i+1);
plot(T,Reference,'r',T,Temp,'b')
xlabel('Time (sec)')
legend('Desired','Simulated')

[Torsten]

Maybe you know why the control does not converge.

k = 20;
Times = linspace(0,0.25,k+1);
dt = 0.25/k;
T_desired = 300.0;      % desired output, or reference point
Kp = 1;                 % proportional term Kp
Ki = 0.5005;            % Integral term Ki
Kd = 0.001;             % derivative term Kd
% pre-assign all the arrays to optimize simulation time
Prop(1:k+1) = 0; Der(1:k+1) = 0; Int(1:k+1) = 0; I(1:k+1) = 0;
PID(1:k+1) = 0;
Temp_avg(1:k+1) = 0;
Temp_avg(1) = 293.15;
voltage (1:k+1) = 0;
Error(1:k+1) = 0;
Error(1) = T_desired - 293.15;
state1(1:k+1) = 0; STATE1(1:k+1) = 0;
state2(1:k+1) = 0; STATE2(1:k+1) = 0;
x1=linspace(0,55,100);
x2=linspace(55.01,59,100);
x3=linspace(59.01,589,100);
x=cat(2,x1,x2,x3); 
ic = @(x) incondfull1(x);
m = 1;
%Times = linspace(0,Time,k+1);
for i=1:k
    C.V=voltage(i);
    eqn = @(x,t,T,dTdx) heateqfull1(x,t,T,dTdx,C);
    sol = pdepe(m,eqn,ic,@boundarycondfull1,x,[Times(i),(Times(i)+Times(i+1))/2,Times(i+1)]);
    Temp_avg(i+1) = 2/(x(110)^2-x(10)^2)*trapz(x(10:110),sol(end,10:110,1).*x(10:110));
    Error(i+1)= T_desired - Temp_avg(i+1);
    Prop(i+1) = Error(i+1);% error of proportional term
    Der(i+1)  = (Error(i+1) - Error(i))/dt; % derivative of the error
    Int(i+1)  = (Error(i+1) + Error(i))*dt/2; % integration of the error
    I(i+1)    = sum(Int); % the sum of the integration of the error    
    PID(i+1)  = Kp*Prop(i) + Ki*I(i+1)+ Kd*Der(i); % the three PID terms
    %% You can replace the follwoing five lines with your system/hardware/model
    STATE1(i+1) = sum(PID); % sum PID term to calculate the first integration
    state2(i+1) = (STATE1(i+1) + STATE1(i))*dt/2; % output after the first integrator
    STATE2(i+1) = sum(state2); % sum output of first integrator to calculate the second integration
    voltage(i+1) = (STATE2(i+1) + STATE2(i))*dt/2; 
end
Reference = T_desired*ones(1,numel(Times));
plot(Times,Reference,'r',Times,Temp_avg,'b');
xlabel('Time (sec)')
legend('Desired','Simulated')

function [c,f,s] = heateqfull1(x,t,T,dTdx,C)
           %%---------- universal parametrs -------------
  tm= 1.7;                            % thickness of membrane in micrometer 
  th= 0.3;                            % thickness of heater in micrometer
  rh= 55;                             % radius of heater in micrometer
  rm= 589;                            % radius of membrane in micrometer 
  wh= 4;                              % width of heater in micrometer 
  r1= 55;                             % radius of region 1 in micrometer 
  r2= 59;                             % radius of region 2 in micrometer 
  r3= 589;                            % radius of region 3 in micrometers
  k= 1.61*10^-5;                      % effective thermal conductivity of Si3N4 and SiO2 in W/µmK 
  ka= 2.623 * 10^-8;                  % thermal conductivity of air in W/µmK 
  kw= 0.598 * 10^-6;                  % thermal conductivity of air in W/µmK 
  hc= 2*10^-10;                       % effective heat trasnfer co-efficent in W/µm2K
  rho= 2.76*10^-12;                   % denisty of membrane g/um3
  shc= 0.68812;                       % specific heat capacity of membrane in J/(g*K)
  Ta= 293.15;                         % room temperature in K 
  kpt= 2.95*10^-5;                    % effective thermal conductivity of platinum in W/µmK

     %%__ resistance of heater depending upon thickness and resistivity ___
         
  rrho= 0.105;                        % resistivity of platinum in ohm um ;  1.05*10*-7 ohm m 2014 paper 
  St=1;                               %2013
  Sc=1;                               %2013
  Wc=24;                              % extrernal wide width contact (nc-ext *wti=12*2=24 )2013
  Sj=20;                              % length of external wide width contact (ns-ext *wti=10*2=20)2013
  Sa=6;                               % 2015 paper supporting documnet as this parameter is not in 2013 paper 
  Lc=503;                             % 2015 paper supporting documnet as this parameter is not in 2013 paper
  th=0.3;                             % thickness of heater in um from 2015 paper 
  Ro=rrho*([(pi*((5*rh)+(6*St)+((25/4)*wh)-(4*Sc)+(2*Sj))/(wh*th))]+[(2*Lc)/(Wc*th)]);
  alpha = 3.927*10^-3;                % temperature coefficient of resistivity in 1/K
  V =C.V;                           % voltage for input in Volts 

         %%--------  specific user defined parameters-----------
  ci= (rho*shc)/k;            
  si=(2*hc)/(k*tm);
  sii=1/(2*kpt*pi*r1*wh*tm*Ro);

  c = ci;
  if x>=0 & x<=55
    s=-si *(T-293.15);
    f = dTdx;
  elseif x>55 & x<=59
    s = -(si*(T-293.15))+ [V^2 *(1-(alpha*(T-293.15)))]/sii ;
    f = dTdx;
  else
    s=-si *(T-293.15);
    f = dTdx;
  end
end

function [pl,ql,pr,qr] = boundarycondfull1(xl,Tl,xr,Tr,t)
  pl = 0;
  ql = 1;
  pr = Tr-293.15  ;
  qr = 0;
end

function T0 = incondfull1(x)
  T0 = 293.15;
end