System Identification เซอร์โวไฮดรอลิก และแบบจำลองทางพลศาสตร์การควบคุมแบบปิด โดย Matlab/Simulink -

เซอร์โวไฮดรอลิก เป็นอุปกรณ์ที่มีใช้อย่างแพร่หลายในอุตสาหกรรมหนัก เพราะด้วยคุณสมบัติของระบบไฮดรอลิกที่ให้กำลังสูง สามารถยกสิ่งของน้ำหนักตั้งแต่ไม่กี่กิโลกรัม ถึงหลายสิบตัน

เซอร์โวไฮดรอลิก ค่อนข้างมีราคาแพงโดยมาแล้วผู้ผลิตจะสร้างสำเร็จรูป เราสามารถนำไปใช้งานได้ทันที แต่ในบล๊อกนี้จะเข้าไปเจาะเบื้องหลังในการสร้างระบบควบคุมของเซอร์โวไฮดรอลิก โดยสามารถเขียนไดอะแกรมได้ดังนี้

จากรูปที่ 1 จะเห็นว่าสามารถแบ่งอุปกรณ์ได้เป็น 2 ส่วน

ส่วน A จะเป็นวาล์วไฮดรอลิก โซลินอยไฟฟ้า รวมถึงอุปกรณ์อื่นๆ ที่ทำหน้าที่ควบคุม ปริมาตรของน้ำมันไหลเข้า-ไหลออก กระบอกสูบ โดยในส่วนนี้เราจะใช้วิธีการ System Identification การประมาณ Transfer function ของระบบที่เราไม่ทราบสมการพลศาสตร์ (System Indentification) ด้วยสมการอนุพันธ์อันดับสอง

ส่วน B จะเป็นชุดกระบอกสูบ ที่ต้องมีกำลังสูงเพื่อทำหน้าที่ ยก กด ดัด และอื่นๆ ที่ต้องใช้กำลังสูง ในส่วนนี้เราจะใช้วิธีหาสมการพลศาสตร์โดยใช้ความรู้ทาง Fluid Mechanics

เราจะเริ่มที่การทำ System Identifation ส่วน A ก่อน

จากการทดสอบเซอร์โววาล์ว พบว่าเมื่อป้อนแรงดันไฟฟ้าควบคุม จาก 0- 10 V จะทำให้ Spool Valve เคลื่อนที่ตามสมการดังนี้

\(y=2.5\times(10^{-5}) E_{in}\)

(1)

ทำการทดสอบป้อนสัญญาณ Step Unit ขนาด 10V ที่ t = 0.5 Sec. และ หยุดป้อนสัญญาณ ที่ t=1.0 Sec. แล้วบันทึกการเคลื่อนที่ของ Spool Valve ด้วยเครื่องมือวัด แสดงได้ดังรูปที่ 2

จากรูปที่ 2 เราพบข้อมูลสำคัญดังนี้

พบว่าวาล์วมีการเคลือนที่จากจุดศูนย์ไปยังตำแหน่ง y=0.25mm
จากผลการทดลองพบว่าค่า Settling Time = 12.7 ms (0.0127 Sec.)
จากผลการทดลองพบว่าระบบแทบไม่เกิด Over shoot ดังนั้น ค่า Dampling Ratio น่าจะค่อนข้างสูง โดยมีค่าระหว่าง 0.8 ถึง 1.0+

Transfer function ของระบบที่มีอนุพันธ์อันดับสอง สามารถแสดงได้ดังบล๊อกไดอะแกรม

G(s) เป็นแบบจำลองที่เราจะทำการ System Identificaion โดยสามารถเขียนได้ดังนี้

\(G(s)=\frac{Y(s)}{U(s)}=\frac{\omega _{n}^{2}}{s^{2}+2\delta \omega _{n}s+\omega _{n}^{2}}\)

(2)

จากความสัมพันธ์ที่ได้จากสมการ (1) จะทำให้เราทราบว่า

\(K_{v}=2.5\times 10^{-5}\: \frac{m}{V}\)

(3)

จากบทความก่อนหน้าเราทราบว่าในระบบพลศาสตร์ที่มีอนุพันธ์อันดับสองจะมีค่า Settling Time ดังนี้

\(T_{s}=\frac{4}{\delta \omega _{n}}\)

(4)

จัดรูปสมการ (4) ทำให้เราทราบว่า

\(\delta \omega _{n}=\frac{4}{T_{s}}=\frac{4}{0.0127}=314.96\)

(5)

ยังคงเหลือค่า Dampling ration \(\delta\) ที่เรายังไม่ทราบค่าแน่ชัด

จากบทความก่อนหน้าเรื่อง การประมาณ Transfer function ของระบบที่เราไม่ทราบสมการพลศาสตร์ (System Indentification) ด้วยสมการอนุพันธ์อันดับสอง

ผู้เขียนได้ทดสอบรันสคริป และเปลี่ยนค่า =0.93 และจากสมการ (5) ทำให้ได้ค่า =338 พบว่าให้ผลตอบสนองใกล้เคียงกับผลการทดลองในรูปที่ 2 มากที่สุด

จากสมการ (2) จะเขียน Transfer function ได้เป็น

\(G(s)=\frac{Y(s)}{U(s)}=\frac{114,244}{s^{2}+2(0.93)(338)s+114,224}\)

(6)

หาแบบจำลองทางพลศาสตร์ของ ส่วน B

ให้ปริมาตรน้ำมันไหลผ่าน Spool valve เข้ากระบอกสูบแสดงได้ด้วยสมการ

\(Q_{1}=C_{d}A_{v}\sqrt{\frac{2}{\rho }(P_{s}-P_{1})}\)

(7)

และในขณะเดียวกันน้ำมันจะไหลออกจากกระบอกสูบ ผ่าน Spool valve กลับลงถัง แสดงได้ด้วยสมการ

\(Q_{2}=C_{d}A_{v}\sqrt{\frac{2}{\rho }(P_{2}-P_{r})}\)

(8)

สมมติให้การไหลของน้ำมันเป็นแบบต่อเนื่อง(steady state) และ น้ำมันไม่สามารถบีบอัดตัวได้ (incompressible) จะได้ว่า \(Q_{1}=Q_{2}\) จะได้สมการ

\(P_{s}-P_{1}=P_{2}-P_{r}\)

(9)

ให้ \(\Delta P=P_{1}-P_{2}\) และ \(\Delta P> 0\) เพราะน้ำมันไหลจาก chamber 2 กลับลงอ่าง สามารถจัด สมการ (9) ใหม่ได้เป็น

\(P_{s}=P_{1}+P_{1}-\Delta P-P_{r}\)

(10)

ปกติแล้ว ความดัน \(P_{r}\) จะน้อยกว่าความดันตัวอื่นๆ มาก จึงสามารถตัดเทอมนี้ทิ้งได้ จัดรูปสมการ (10) ได้เป็น

\(P_{1}=\frac{P_{s}+\Delta P}{2}\)

(11)

แทนสมการ (11) กลับลงใน สมการ (7) จะได้

\(Q_{1}=C_{d}A_{v}\sqrt{\frac{2}{\rho }\left ( P_{s}-\frac{P_{s}+\Delta P}{2} \right )}\)

(12)

จากความสัมพันธ์พื้นที่หน้าตัด \(A_{v}=hy \) สมการ (12) เขียนใหม่ได้เป็น

\(Q_{1}=C_{d}hy\sqrt{\frac{2}{\rho }\left ( P_{s}-\frac{P_{s}+\Delta P}{2} \right )}\)

(13)

ทำการ Linearized สมการ (13) รอบๆ จุดทำงาน \(y^{\ast }\) และ \(\Delta P^{\ast }\)

\(\delta Q_{1}=\frac{\partial f}{\partial y}\delta y+\frac{\partial f}{\partial \Delta P}\delta (\Delta P)\)

(14)

จะได้อนุพันธ์ย่อยของสมการ (14) ดังนี้

\(\frac{\partial f}{\partial y}=C_{d}h\sqrt{\frac{P_{s}-\Delta P}{\rho }}\)

(15)

และ

\(\frac{\partial f}{\partial \Delta P}=-\frac{C_{d}hy}{2\varrho }\left ( \frac{P_{s}-\Delta P}{\varrho } \right )^{-\frac{1}{2}}\)

(16)

ณ จุดที่ทำการ linearized ให้เป็นสภาวะที่ยังไม่เกิดการเปลี่ยนแปลง nominal state \(y^{\ast }\)=0 และ \(\Delta P^{\ast }\)=0 (no flow) จะได้

\(\frac{\partial f}{\partial y}=C_{d}h\sqrt{\frac{P_{s}}{\varrho }}\)

(17)

และ

\(\frac{\partial f}{\partial \Delta P}=0\)

(18)

สมการ (14) ที่ทำการแปลงให้เป็นเชิงเส้นแล้วจะสามารถเขียนใหม่ได้เป็น

\(\delta Q_{1}=C_{d}h\sqrt{\frac{P_{s}}{\varrho }}\delta y\)

(19)

จากสภาวะรอบๆ จุด ที่ทำการแปลงเป็นเชิงเส้น \(\delta Q_{1}=Q_{1}-Q_{1}^{\ast }\) และ \(\delta y=y-y^{\ast }\) รวมถึงสมมติฐานที่ให้สภาวะ ณ จุดนี้ยังไม่เกิดการเปลี่ยนแปลง \(Q_{1}^{\ast }\)=0 , \(y^{\ast }\)=0 เขียนสมการ (19) ใหม่ได้เป็น

\(Q_{1}=C_{d}h\sqrt{\frac{P_{s}}{\varrho }}y\)

(20)

อนุพันธ์ของปริมาตร \(\dot{V_{1}}=A\dot{x}\) สมการ (20) จะกลายเป็น

\(C_{d}h\sqrt{\frac{P_{s}}{\varrho }}y=A\dot{x}\)

(21)

จัดรูปสมากร (21) เราจะได้สมการพลศาสตร์ของส่วน B

\(\dot{x}=K_{HA}\: y\)

(22)

โดยที่ \(K_{HA}=\frac{C_{d}h}{A}\sqrt{\frac{P_{s}}{\varrho }}\)

ทำการอินทริเกรท สมการ (22) จะได้

\(x(t)=x_{0}+K_{HA}\int ydt\)

(23)

โดยที่ \(x_{0}\) คือตำแหน่งเริ่มต้นของกระบอกสูบ

ทำการจำลงอการทำงานใน Matlab/Simulink โดยการสร้างคอนโทรลบล๊อกไดอะแกรมดังภาพ

และทำการรับ Script ดังนี้

%%  PARAMETERS
m = 12; % kg Piston and load mass
b = 215; % Viscous friction coefficient, b, N.s/m
A = 6.78*(10^(-4)); % Piston area, A,m2
V0 = 1.58*(10^(-4)); % Minimum chamber volume, V0,m3
L = 0.68; % Piston stroke, L, m
P_s = 18*10^6; % Supply pressure, PS, MPa
P_r = 0.12133*(10^(6)); %Reservoir (drain) pressure, Pr, MPa
beta = 729*(10^(6)); % Fluid bulk modulus, ?  ,MPa
Cd = 0.72; % Valve discharge coefficient, Cd
rho = 899; % Fluid density, ?  , kg/m3
h = 0.0094; % Valve orifice height, h, m
wn = 338; % Solenoid–valve undamped natural frequency, ?n , rad/s
zeta = 0.93; % Solenoid–valve damping ratio, ? 
Kv = 25*(10^(-6)); % Solenoid–valve DC gain, Kv  mm/V
K_HA = ((Cd*h)/A)*sqrt(P_s/rho); % hydraulic actuator gain
KP=1550;
K_LF=2000;
%%
sysG=tf(4034,[1 628 114244 0 ]);
rlocus(sysG)
title('Root Locus')
%%
Kp=1550
Close_Loop_poles=rlocus(sysG,Kp)
%%
open("Proportional.slx")
sim("Proportional.slx")

จะได้ผลลัพธ์ดังนี้