現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
Pdepe solver be run in loop for designing a control scheme.
2 ビュー (過去 30 日間)
古いコメントを表示
maria atiq
2024 年 1 月 30 日
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 件のコメント
maria atiq
2024 年 1 月 30 日
When i try to give index to the sol() term . Matlab gives the error that it is illegal to index this function.
Torsten
2024 年 1 月 30 日
編集済み: Torsten
2024 年 1 月 30 日
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
2024 年 1 月 30 日
編集済み: maria atiq
2024 年 1 月 30 日
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
2024 年 1 月 30 日
I think we are unable to help you without some code that clarifies your problem.
maria atiq
2024 年 1 月 30 日
編集済み: maria atiq
2024 年 1 月 30 日
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
2024 年 1 月 30 日
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
2024 年 1 月 30 日
What does the index i stand for ? The average solution for T over all x-values at time (i*dt) ?
maria atiq
2024 年 1 月 30 日
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
2024 年 1 月 30 日
編集済み: Torsten
2024 年 1 月 30 日
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
2024 年 1 月 30 日
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
2024 年 2 月 7 日
編集済み: maria atiq
2024 年 2 月 7 日
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
2024 年 2 月 7 日
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
2024 年 2 月 7 日
編集済み: Torsten
2024 年 2 月 7 日
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
2024 年 2 月 7 日
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
2024 年 2 月 7 日
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
2024 年 2 月 7 日
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
2024 年 2 月 7 日
編集済み: maria atiq
2024 年 2 月 7 日
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
2024 年 2 月 7 日
編集済み: Torsten
2024 年 2 月 7 日
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
2024 年 2 月 8 日
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
2024 年 2 月 8 日
編集済み: Torsten
2024 年 2 月 8 日
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
2024 年 2 月 10 日
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
2024 年 2 月 10 日
It would mean alot if this issue can be sorted as i have been stuck at it for way too long. Kind regards.
Torsten
2024 年 2 月 10 日
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
2024 年 2 月 10 日
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
2024 年 2 月 10 日
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
2024 年 2 月 13 日
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
2024 年 2 月 14 日
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
回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Heat and Mass Transfer についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の MathWorks のサイトは、お客様の地域からのアクセスが最適化されていません。
南北アメリカ
- América Latina (Español)
- Canada (English)
- United States (English)
ヨーロッパ
- 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)
アジア太平洋地域
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)