From 8fca4cd8c28d990c0ff4009159804002fa95e6f9 Mon Sep 17 00:00:00 2001 From: Jens Ahrensfeld Date: Mon, 10 Sep 2018 16:51:37 +0000 Subject: [PATCH] - updated git-svn-id: http://moon:8086/svn/projects/HendiControl@5 fda53097-d464-4ada-af97-ba876c37ca34 --- matlab/brew_eval.m | 184 +++++++++++++++++++++++++++++++++------------ 1 file changed, 134 insertions(+), 50 deletions(-) diff --git a/matlab/brew_eval.m b/matlab/brew_eval.m index d2e35da..e3206a9 100755 --- a/matlab/brew_eval.m +++ b/matlab/brew_eval.m @@ -18,14 +18,13 @@ function brew_eval(varargin) % -% brew_eval('M', 20.0, 'cont_ctrl', 0, 'kp', 0.1, 'ki', 0.1, 'kd', 100, 'T', 1800, 'mass_kleak', 0.2, 'P_q', 100, 'mass_delay', 100, [50 60 70 80]); +% brew_eval('M', 20.0, 'cont_ctrl', 1, 'kp', 2.0, 'ki', 0.02, 'kd', 400, 'T', 1800, 'mass_kleak', 0.5, 'P_q', 100, 'mass_delay', 60, 'dt', 1.00, [30 50 70]); -WITH_LOWPASS = 0; -WITH_RUNNING_MEAN = 0; WITH_RATE_CONTROLLER = 1; +WITH_HYST = 0.0; -params = struct ('cont_ctrl',0, 'C',4.19e3, 'M',10, 'mass_kleak',0.5, 'mass_delay',50, 'P_min',0, 'P_max',3000, 'P_q',100, 'T',100, 'dt',1.0, 'Td',10.0, 'theta_lock',0.5, 'theta_amb',20.0, 'kp',10, 'ki',0.04, 'kd',500, 'pid_kleak',0.1, 'k_noise', 0.0, 'k_noise2', 0.0, 'f_noise2', 0.01); -units = struct ('cont_ctrl', '', 'C', 'J/(kg*K)', 'M', 'kg', 'mass_kleak', '', 'mass_delay', 's', 'P_min', 'W', 'P_max', 'W', 'P_q','W', 'T', 's', 'dt', 's', 'Td', 's', 'theta_lock', '°C', 'theta_amb', '°C', 'kp', '', 'ki', '', 'kd', '', 'pid_kleak', '', 'k_noise', '', 'k_noise2', '', 'f_noise2', 'Hz'); +params = struct ('cont_ctrl',0, 'C',4.19e3, 'M',10, 'mass_kleak',0.5, 'mass_delay',50, 'P_min', 0, 'P_max',3000, 'P_q',100, 'T',100, 'dt',1.0, 'Td',10.0, 'theta_lock', 0.5, 'theta_amb',20.0, 'kp',10, 'ki',0.04, 'kd',500, 'pid_kleak',0.001, 'k_noise', 0.0, 'heatrate', 1.0); +units = struct ('cont_ctrl', '', 'C', 'J/(kg*K)', 'M', 'kg', 'mass_kleak', '', 'mass_delay', 's', 'P_min', 'W', 'P_max', 'W', 'P_q','W', 'T', 's', 'dt', 's', 'Td', 's', 'theta_lock', '°C', 'theta_amb', '°C', 'kp', '', 'ki', '', 'kd', '', 'pid_kleak', '', 'k_noise', '', 'heatrate', '°C/min'); % Parse parameters names = fieldnames(params); @@ -48,45 +47,121 @@ theta_last = params.theta_amb; thetas = varargin{nargin} heatRate = 0; heatRate_f = 0; -heatRateMeasureInterval = 10.0; % seconds +heatRateMeasureInterval = 1.0; % seconds heatRateMeasureCount = 0; % seconds n = 0; nSteps = params.P_max / params.P_q; N = 2*fix(length(thetas)*params.T/params.dt); -maz = 17; -lp_N = 17; -[lp_b, lp_a] = butter (lp_N-1, 0.125); +maz = 7; +lp_N = 5; +[lp_b, lp_a] = butter (lp_N-1, 0.25); lp_z = zeros(1, length(lp_b)-1); lp_rate_N = 17; [lp_rate_b, lp_rate_a] = butter (lp_rate_N-1, 0.125); lp_rate_z = zeros(1, length(lp_rate_b)-1); +P_cont = 0; +y_rate = 0; +y = 0; + +stbl_rate_z = []; +stbl_hold_z = []; + +ovenState = 'HOLD'; +controllerState = 'TRK'; + +curr = struct('kp', 0, 'ki', 0, 'kd', 0); +curr_rate = struct('kp', 0, 'ki', 0, 'kd', 0); + +err0_hold = 0; +err0_rate = 0; +a_err = 0.5; +P_max = 0; +heatRate_soll = 0; for k=1:length(thetas) theta_soll = thetas(k); hold_theta = 0; - temp_ok = 0; + P_max = 0; + printf('%d: Target temperatur %d °C\n', k, theta_soll); while n < N, - err = theta_ist - theta_soll; - [y, yc, s_pid] = jpid(s_pid, params.kp, params.ki*params.dt, params.kd*params.dt, 1-params.pid_kleak, -err); - s_pid.y_max = 50; - - P_rate = params.P_max; - % Rate controller - err_rate = heatRate - 1.0; - [y_rate, yc_rate, s_pid_rate] = jpid(s_pid_rate, 0.4*params.dt, 0.001*params.dt, 8.0*params.dt, 0.999, -err_rate); - s_pid_rate.y_max = 4000; - if WITH_RATE_CONTROLLER - if (hold_theta == 0) - P_rate = min(params.P_max, params.P_max*y_rate); - end + err1 = theta_ist - theta_soll; + err1_rate = heatRate - heatRate_soll; + err0_hold = (1-a_err)*err0_hold + a_err*err1; + err0_rate = (1-a_err)*err0_rate + a_err*err1_rate; + [err_rate_stbl, stbl_rate_z] = isStable(stbl_rate_z, err0_rate, 0.05); + [err_hold_stbl, stbl_hold_z] = isStable(stbl_hold_z, err0_hold, 0.05); + ovenState_n = ovenState; + controllerState_n = controllerState; + switch ovenState + case {'HEAT'} + heatRate_soll = params.heatrate; + curr_rate.kp = 4.0; + curr_rate.ki = 0.02; + curr_rate.kd = 4.0; + P_max = 0.5*P_max + 0.5*P_cont; + switch controllerState + case {'ACQ'} + if err_rate_stbl && (abs(err0_rate) < 0.1) + controllerState_n = 'TRK'; + end + case {'TRK'} + curr_rate.kp = 8.0; + curr_rate.ki = 0.03; + curr_rate.kd = 8.0; + if abs(err0_rate) > 0.5 + controllerState_n = 'ACQ'; + end + endswitch + if abs(err0_hold) < 1.0 + ovenState_n = 'HOLD'; + controllerState_n = 'ACQ'; + end + case {'HOLD'} + heatRate_soll = 0; + curr.kp = params.kp; + curr.ki = params.ki; + curr.kd = params.kd; + switch controllerState + case {'ACQ'} + if err_hold_stbl && (abs(err0_hold) < 0.1) + controllerState_n = 'TRK'; + end + case {'TRK'} + curr.kp = params.kp; + curr.ki = params.ki; + curr.kd = params.kd; + if abs(err0_hold) > 0.2 + controllerState_n = 'ACQ'; + end + endswitch + if err0_hold < -1.0 + if WITH_RATE_CONTROLLER + ovenState_n = 'HEAT'; + end + controllerState_n = 'ACQ'; + end + endswitch + if strcmp(controllerState, controllerState_n) == 0 || strcmp(ovenState, ovenState_n) == 0 + printf('%d: State %s::%s -> %s::%s\n', n, ovenState, controllerState, ovenState_n, controllerState_n); end - % Rate controller + ovenState = ovenState_n; + controllerState = controllerState_n; - P_cont = max(params.P_min, min(P_rate, P_rate*y)); - P_cont_q = P_cont; %fix(P_norm * nSteps)*params.P_q; + % Hold controller + if strcmp(ovenState, 'HEAT') + [y, yc, s_pid] = jpid(s_pid, curr_rate.kp, curr_rate.ki*params.dt, curr_rate.kd/params.dt, 1-params.pid_kleak, -err0_rate); + P_cont = max(params.P_min, min(params.P_max, params.P_max*y)); + s_pid.y_max = 80; + else + [y, yc, s_pid] = jpid(s_pid, curr.kp, curr.ki*params.dt, curr.kd/params.dt, 1-params.pid_kleak, -err0_hold); + P_cont = max(params.P_min, min(P_max, params.P_max*y)); + s_pid.y_max = 80; + end + + P_cont_q = fix(P_cont/params.P_q + 0.5)*params.P_q; if params.cont_ctrl P = P_cont_q; else @@ -95,43 +170,34 @@ for k=1:length(thetas) P = params.P_max * heat; end [theta, s_mass] = mass(s_mass, params.dt, params.C, params.M, params.mass_kleak, P, params.mass_delay, theta_ist-params.theta_amb); - theta_ist = params.theta_amb + theta + params.k_noise*randn()/sqrt(12) + params.k_noise2*cos(params.f_noise2*2*pi*n*params.dt); + theta_ist = params.theta_amb + theta + params.k_noise*randn()/sqrt(12); - if WITH_RUNNING_MEAN - [theta_ist, maz] = movingavg(theta_ist, maz); - end - if WITH_LOWPASS - [theta_ist, lp_z] = filter(lp_b, lp_a, theta_ist, lp_z); - end if heatRateMeasureCount <= 0 heatRate = 60*(theta_ist-theta_last)/heatRateMeasureInterval; heatRateMeasureCount = heatRateMeasureInterval; - theta_last = theta_ist ; - else - heatRateMeasureCount = heatRateMeasureCount - params.dt; + theta_last = theta_ist; end + heatRateMeasureCount = heatRateMeasureCount - params.dt; [heatRate_f, lp_rate_z] = filter(lp_rate_b, lp_rate_a, heatRate, lp_rate_z); n = n + 1; theta_ist = theta_ist; heatRate_(n) = heatRate_f; theta_(n) = theta_ist; - err_(n) = err; + err_(n) = err0_hold; p_(n) = P; pc_(n) = P_cont; pcq_(n) = P_cont_q; y_(n) = y; - P_rate_(n) = P_rate; - err_rate_(n) = err_rate; + y_rate_(n) = y_rate; + err_rate_(n) = err0_rate; + if hold_theta > 0 hold_theta = hold_theta - 1; if hold_theta == 0 - disp('Temperature next') break; end - elseif abs(err) < params.theta_lock + elseif abs(err0_hold) < params.theta_lock hold_theta = params.T/params.dt; - temp_ok = 1; - disp('Temperature is OK!') end end end @@ -141,17 +207,17 @@ t = (0:n-1)*params.dt; subplot(3,1,1) plot(t, theta_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Temperature'); subplot(3,1,2) -plot(t, p_, 'b-', t, y_, 'g-', pc_, 'r-'); grid; xlabel('t/s'); ylabel('P/W'); title('Power Control'); legend('P','y','P_{cont}'); +plot(t, p_, 'b-', t, pcq_, 'g-', t, pc_, 'r-'); grid; xlabel('t/s'); ylabel('P/W'); title('Power Control'); legend('P','y','P_{cont}'); subplot(3,1,3) -plot(t, err_); grid; xlabel('t/s'); ylabel('err/K'); title('Error'); legend('Error'); +plot(t, err_, t, err_rate_); grid; xlabel('t/s'); ylabel('err/K'); title('Error'); legend('Error_{Hold}', 'Error_{Rate}'); figure; subplot(3,1,1) -plot(t, heatRate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Aufheizrate'); +plot(t, heatRate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Heat_{rate}'); subplot(3,1,2) -plot(t, P_rate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('P_rate'); +plot(t, y_rate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('y_{rate}'); subplot(3,1,3) -plot(t, err_rate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Error'); +plot(t, err_rate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Error_{rate}'); function [theta, s] = mass(si, dt, C, M, L, P, Td, dTheta) s = si; @@ -160,10 +226,10 @@ if ~isstruct(si) if Td > 0 alpha = dt/Td; end - s = struct('e', 0, 'a', alpha, 'x', 0, 'gain', 0.7); + s = struct('e', 0, 'a', alpha, 'x', 0, 'gain', 0.8); end -s.e = s.e * (1-((L*dTheta)*dt)/(M*C)); +s.e = s.e * (1-((L*(0+dTheta))*dt)/(M*C)); s.x = (1-s.a)*s.x + s.gain*s.a*P*dt; s.e = s.e + s.x; theta = s.e/(M*C); @@ -194,6 +260,24 @@ end endfunction +function [y, z] = isStable(zi, v, maxErr) +z = zi; +if ~isstruct(zi) + z = struct('v', 0, 'cnt', 3, 'cnt_reload', 3); +end + + +if z.cnt > 0 + z.cnt = z.cnt - 1; +endif + +if abs(v - z.v) <= maxErr + z.cnt = z.cnt_reload; +endif + +y = (z.cnt == 0); +endfunction + function print_parameters(params) align = 32; names = fieldnames(params);