You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Pdepe solver be run in loop for designing a control scheme.
2 views (last 30 days)
Show older comments
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.
35 Comments
maria atiq
on 30 Jan 2024
When i try to give index to the sol() term . Matlab gives the error that it is illegal to index this function.
Torsten
on 30 Jan 2024
Edited: Torsten
on 30 Jan 2024
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
on 30 Jan 2024
Edited: maria atiq
on 30 Jan 2024
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
on 30 Jan 2024
I think we are unable to help you without some code that clarifies your problem.
maria atiq
on 30 Jan 2024
Edited: maria atiq
on 30 Jan 2024
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
on 30 Jan 2024
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
on 30 Jan 2024
What does the index i stand for ? The average solution for T over all x-values at time (i*dt) ?
maria atiq
on 30 Jan 2024
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);
Torsten
on 30 Jan 2024
Edited: Torsten
on 30 Jan 2024
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
on 30 Jan 2024
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
maria atiq
on 7 Feb 2024
Edited: maria atiq
on 7 Feb 2024
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
on 7 Feb 2024
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
maria atiq
on 7 Feb 2024
Edited: Torsten
on 7 Feb 2024
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')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1610611/image.png)
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
Torsten
on 7 Feb 2024
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
on 7 Feb 2024
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.
Torsten
on 7 Feb 2024
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
on 7 Feb 2024
Edited: maria atiq
on 7 Feb 2024
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
on 7 Feb 2024
Edited: Torsten
on 7 Feb 2024
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
on 8 Feb 2024
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
on 8 Feb 2024
Edited: Torsten
on 8 Feb 2024
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
on 10 Feb 2024
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
on 10 Feb 2024
It would mean alot if this issue can be sorted as i have been stuck at it for way too long. Kind regards.
Torsten
on 10 Feb 2024
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
on 10 Feb 2024
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
on 10 Feb 2024
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.
maria atiq
on 13 Feb 2024
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')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1615773/image.png)
Torsten
on 14 Feb 2024
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')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1617138/image.png)
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
Answers (0)
See Also
Categories
Find more on Heat and Mass Transfer in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)