Optimalisasi sistem kontrol LQR

Pendahuluan

Beberapa artikel [1,2,3] telah dipublikasikan di Habr yang secara langsung atau tidak langsung berkaitan dengan topik ini. Dalam hal ini, orang tidak dapat gagal untuk mencatat publikasi [1] dengan judul "Matematika pada Jari: Regulator Linier-Kuadratik", yang populer menjelaskan prinsip operasi pengontrol LQR optimal.

Saya ingin melanjutkan topik ini, setelah memeriksa aplikasi praktis dari metode optimasi dinamis, tetapi sudah pada contoh nyata menggunakan Python. Pertama, beberapa kata tentang terminologi dan metode optimasi dinamis.

Metode optimasi dibagi menjadi statis dan dinamis. Objek kontrol berada dalam keadaan gerakan terus menerus di bawah pengaruh berbagai faktor eksternal dan internal. Oleh karena itu, evaluasi hasil kontrol diberikan untuk waktu kontrol T, dan ini adalah tugas optimasi dinamis.

Menggunakan metode optimisasi dinamis, masalah yang terkait dengan distribusi sumber daya terbatas selama periode waktu tertentu diselesaikan, dan fungsi tujuan ditulis sebagai fungsional integral.

Alat matematika untuk memecahkan masalah tersebut adalah metode variasional: kalkulus variasi klasik, prinsip L.S. maksimum. Pontryagin dan pemrograman dinamis R. Bellman.

Analisis dan sintesis sistem kontrol dilakukan dalam ruang waktu, frekuensi, dan keadaan. Analisis dan sintesis sistem kontrol di ruang negara dimasukkan ke dalam kurikulum, namun, teknik yang disajikan dalam materi pelatihan menggunakan pengontrol SQR dirancang untuk menggunakan Matlab dan tidak mengandung contoh praktis analisis.

Tujuan dari publikasi ini adalah untuk mempertimbangkan metode analisis dan sintesis sistem kontrol linier di ruang negara menggunakan contoh terkenal mengoptimalkan sistem kontrol drive listrik dan pesawat terbang menggunakan bahasa pemrograman Python.

Pembuktian matematis dari metode optimasi dinamis

Untuk optimasi, kriteria amplitudo (MO), simetris (CO) dan kompromi optimum (KO) digunakan.

Ketika memecahkan masalah optimisasi dalam ruang keadaan, sistem stasioner linier diberikan oleh persamaan vektor-matriks:

$$$\ dot {x} = \ frac {dx} {dt} = A \ cdot x + B \ cdot u; y = C \ cdot x,$ (1)

Kriteria integral untuk konsumsi energi kontrol minimum dan kecepatan maksimum ditetapkan oleh fungsi:

$$$J_ {x} = \ int_ {0} ^ {\ infty} (x ^ {T} Qx + u ^ {T} Ru + 2x ^ {T} Nu) dt \ rightarrow min,$ (2)

$$$J_ {u} = \ int_ {0} ^ {\ infty} (y ^ {T} Qy + u ^ {T} Ru + 2y ^ {T} Nu) dt \ rightarrow min.$ (3)

Hukum kontrol u adalah dalam bentuk umpan balik linier pada variabel status x atau pada variabel output y:

$$$u = \ pm Kx, u = \ pm Ky$

Kontrol tersebut meminimalkan kriteria kualitas kuadrat (2), (3). Dalam hubungan (1) ÷ (3), Q dan R adalah matriks pasti definitif positif simetris dari dimensi [n × n] dan [m × m], masing-masing; K adalah matriks koefisien konstan dimensi [m × n], yang nilainya tidak terbatas. Jika parameter input N dihilangkan, ini diambil menjadi nol.

Solusi untuk masalah ini, yang disebut masalah linear integral kuadrat optimasi (LQR-optimasi), dalam ruang keadaan ditentukan oleh ekspresi:

$$$u = R ^ {- 1} B ^ {T} Px$

di mana matriks P harus memenuhi persamaan Ricatti:

$$$A ^ {T} P + PA + PBR ^ {- 1} B ^ {T} P + Q = 0$

Kriteria (2), (3) juga bertentangan, karena penerapan istilah pertama membutuhkan sumber daya yang sangat tinggi, dan untuk yang kedua, sumber daya yang sangat rendah. Ini dapat dijelaskan dengan ungkapan berikut:

$$$J_ {x} = \ int_ {0} ^ {\ infty} x ^ {T} Qxdt$ ,

yang merupakan norma $$$Mod (x)$ vektor x, mis. ukuran osilasi dalam proses regulasi, dan, oleh karena itu, mengambil nilai yang lebih kecil dalam transien cepat dengan osilasi lebih sedikit, dan ungkapan:

$$$J_ {u} = \ int_ {0} ^ {\ infty} u ^ {T} Rudt$

adalah ukuran dari jumlah energi yang digunakan untuk kontrol, itu adalah penalti untuk biaya energi sistem.

Matriks bobot Q, R, dan N menentukan batasan koordinat yang sesuai. Jika ada elemen dari matriks ini sama dengan nol, maka koordinat yang sesuai tidak memiliki batasan. Dalam praktiknya, pemilihan nilai matriks Q, R, N dilakukan dengan metode estimasi ahli, uji coba, kesalahan dan tergantung pada pengalaman dan pengetahuan insinyur desain.

Untuk mengatasi masalah tersebut, operator MATLAB berikut digunakan:
$$$\ begin {bmatrix} R, S, E \ end {bmatrix} = lqr (A, B, Q, R, N)$ dan $$$\ begin {bmatrix} R, S, E \ end {bmatrix} = lqry (P_ {s}, Q, R, N)$ yang meminimalkan fungsional (2), (3) oleh vektor negara x atau oleh vektor keluaran y.

Model Obyek Manajemen $$$P_ {s} = ss (A, B, C, D).$ Hasil perhitungan adalah matriks K dari umpan balik optimal sehubungan dengan variabel keadaan x, solusi dari persamaan Ricatti S dan nilai eigen E = eiq (A-BK) dari sistem kontrol loop tertutup.

Komponen fungsional:

$$$J_ {x} = x0 ^ {T} P_ {1} x0, J_ {u} = x0 ^ {T} P_ {2} x0$

di mana x0 adalah vektor kondisi awal; $$$P_ {1}$ dan $$$P_ {2}$ - Matriks yang tidak diketahui yang merupakan solusi dari persamaan matriks Lyapunov. Mereka ditemukan menggunakan fungsi; $$$P_ {1} = lyap (NN, Q)$ dan $$$P_ {2} = lyap (NN, K ^ {T} RK)$ , $$$NN = (A + BK) ^ {T}$

Fitur implementasi metode optimasi dinamis menggunakan Python

Meskipun Python Control Systems Library [4] memiliki fungsi: control.lqr, control.lyap, penggunaan control.lqr hanya dimungkinkan jika modul slycot diinstal, yang merupakan masalah [5]. Ketika menggunakan fungsi lyap dalam konteks tugas, hasil optimasi dalam control.exception.ControlArgument: Q harus menjadi kesalahan matriks simetris [6].

Oleh karena itu, untuk mengimplementasikan fungsi lqr () dan lyap (), saya menggunakan scipy.linalg:

from numpy import* from scipy.linalg import* def lqr(A,B,Q,R): S =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*S)) E= eig(AB*K)[0] return K, S, E P1=solve_lyapunov(NN,Ct*Q*C) 

Omong-omong, Anda seharusnya tidak sepenuhnya meninggalkan kontrol, karena fungsi step (), pole (), ss (), tf (), feedback (), acker (), place () dan lainnya berfungsi dengan baik.

Contoh optimasi LQR dalam drive listrik

(contoh diambil dari publikasi [7])

Pertimbangkan sintesis pengontrol linear-kuadratik yang memenuhi kriteria (2) untuk objek kontrol yang didefinisikan dalam ruang keadaan oleh matriks:

$$

$$A = \ begin {bmatrix} & -100 & 0 & 0 \\ & 143.678 & -16.667 & -195.402 \\ & 0 & 1.046 & 0 \ end {bmatrix}; B = \ begin {bmatrix} 2300 \\ 0 \\ 0 end {bmatrix}; C = \ begin {bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \ end {bmatrix}; D = 0$$

Berikut ini dianggap sebagai variabel keadaan: $$$x_ {1}$ - tegangan konverter, V; $$$x_ {2}$ - Arus motor, A; $$$x_ {3}$ - kecepatan sudut $$$c ^ {- 1}$ . Ini adalah sistem TP-DPT dengan HB: engine R nom = 30 kW; Unom = 220 V; Inom = 147 A ;; $$$\ omega ω$ 0 = 169 $$$c ^ {- 1}$ ; $$$\ omega ω$ maks = 187 $$$c ^ {- 1}$ ; momen resistansi nominal Mnom = 150 N * m; banyaknya arus mulai = 2; thyristor converter: Unom = 230 V; Uy = 10 B; Inom = 300 A; multiplisitas arus lebih pendek = 1.2.

Saat memecahkan masalah, kami menerima matriks Q diagonal. Sebagai hasil pemodelan, ditemukan bahwa nilai minimum elemen matriks R = 84, dan $$$Q = [[0,01,0,0], [0, 0,01,0], [0,0,0,01]].$ Dalam hal ini, proses transisi monoton dari kecepatan sudut mesin diamati (Gbr. 1). Pada R = 840 $$$Q = [[0,01,0,0], [0, 0,01,0], [0,0,0,01]]$ proses transien (Gbr. 2) memenuhi kriteria MO. Perhitungan matriks P1 dan P2 dilakukan pada x0 = [220 147 162].



Fig. 1



Fig. 2
Dengan pemilihan matriks R dan Q yang tepat, dimungkinkan untuk mengurangi arus awal motor ke nilai yang dapat diterima sama dengan (2–2,5) In (Gbr. 3). Misalnya, dengan R = 840 dan
Q = ([[[0,01,0,0]], [0,0,88,0], [0,0,0,01]], nilainya 292 A, dan waktu proses transisi dalam kondisi ini adalah 1,57 dtk.



Gbr.3

Dalam semua kasus yang dipertimbangkan, umpan balik pada tegangan dan arus negatif, dan pada kecepatan, positif, yang tidak diinginkan sesuai dengan persyaratan stabilitas. Selain itu, sistem yang disintesis adalah astatic untuk tugas dan statis untuk beban. Oleh karena itu, kami mempertimbangkan sintesis pengontrol PI di ruang keadaan dengan variabel status tambahan $$$x_ {4}$ - Koefisien transfer integrator.

Kami mewakili informasi awal dalam bentuk matriks:

$$$A = \ begin {bmatrix} & -100 & 0 & 0 & 0 \\ & 143.678 & -16.667 & -195.402 & 0 \\ & 0 & 1.046 & 0 & 0 \\ & 0 & 0 & 1 & 0 \ end {bmatrix}; B = \ begin {bmatrix} 2300 \\ 0 \\ 0 \\ 0 \\ end {bmatrix}; C = eye (4)); D = 0$

Transien tugas yang sesuai dengan kriteria MO diperoleh pada R = 100, q11 = q22 = q33 = 0,001 dan q44 = 200. Gambar 4 menunjukkan transien variabel keadaan yang mengkonfirmasi astatisme sistem dengan tugas dan beban.



Fig. 4
Untuk menentukan matriks K, Pustaka Sistem Kontrol Python memiliki dua fungsi K = acker (A, B, s) dan K = place (A, B, s), di mana s adalah vektor - deretan kutub yang diinginkan dari fungsi transfer sistem kontrol loop tertutup. Perintah pertama hanya dapat digunakan untuk sistem dengan satu input di u untuk n≤5. Yang kedua tidak memiliki batasan seperti itu, tetapi multiplisitas kutub tidak boleh melebihi pangkat matriks B. Contoh penggunaan operator acker () ditunjukkan dalam daftar dan dalam (Gbr. 5).




Fig. 5

Kesimpulan

Implementasi LQR-optimasi drive listrik yang diberikan dalam publikasi, dengan mempertimbangkan penggunaan fungsi lqr () dan lyap () baru, akan memungkinkan kita untuk meninggalkan penggunaan program berlisensi MATLAB ketika mempelajari bagian yang sesuai dari teori kontrol.

Contoh optimasi LQR di pesawat terbang (LA)

(contoh diambil dari publikasi oleh Astrom dan Mruray, Bab 5 [8] dan diselesaikan oleh penulis artikel dalam hal mengganti fungsi lqr () dan mengurangi terminologi menjadi yang diterima secara umum)

Bagian teoretis, teori singkat, optimasi LQR telah dipertimbangkan di atas, jadi mari kita beralih ke analisis kode dengan komentar yang sesuai:

Pustaka yang diperlukan dan fungsi pengontrol LQR.

 from scipy.linalg import* # scipy   lqr from numpy import * # numpy    from matplotlib.pyplot import * #     from control.matlab import * #  control..ss()  control..step() #     lqr() m = 4; #    . J = 0.0475; #        ***2 #́-       #  . r = 0.25; #     . g = 9.8; #     /**2 c = 0.05; #   () def lqr(A,B,Q,R): X =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*X)) E= eig(AB*K)[0] return X, K, E 

Kondisi awal dan matriks dasar model A, B, C, D.

 xe = [0, 0, 0, 0, 0, 0];#   x   () ue = [0, m*g]; #   #   A = matrix( [[ 0, 0, 0, 1, 0, 0], [ 0, 0, 0, 0, 1, 0], [ 0, 0, 0, 0, 0, 1], [ 0, 0, (-ue[0]*sin(xe[2]) - ue[1]*cos(xe[2]))/m, -c/m, 0, 0], [ 0, 0, (ue[0]*cos(xe[2]) - ue[1]*sin(xe[2]))/m, 0, -c/m, 0], [ 0, 0, 0, 0, 0, 0 ]]) #   B = matrix( [[0, 0], [0, 0], [0, 0], [cos(xe[2])/m, -sin(xe[2])/m], [sin(xe[2])/m, cos(xe[2])/m], [r/J, 0]]) #   C = matrix([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]]) D = matrix([[0, 0], [0, 0]]) 

Kami membangun input dan output yang sesuai dengan perubahan langkah demi langkah pada posisi xy. Vektor xd dan yd sesuai dengan kondisi stabil yang diinginkan dalam sistem. Matriks Cx dan Cy adalah output yang sesuai dari model. Untuk menggambarkan dinamika sistem, kami menggunakan sistem persamaan vektor-matriks:

$$

$$\ left \ {\ begin {matrix} xdot = Kapak + Bu => xdot = (A-BK) x + xd \\ u = -K (x-xd) \\ y = Cx \\ \ end {matrix} \ benar$$

Dinamika loop tertutup dapat dimodelkan menggunakan fungsi step (), untuk K \ cdot xd dan K \ cdot xd sebagai vektor input, di mana:

 xd = matrix([[1], [0], [0], [0], [0], [0]]); yd = matrix([[0], [1], [0], [0], [0], [0]]); 

Pustaka 0.8.1 kontrol saat ini hanya mendukung sebagian fungsi SISO Matlab untuk membuat pengontrol umpan balik, oleh karena itu, untuk membaca data dari pengontrol, oleh karena itu diperlukan untuk membuat vektor indeks lat (), alt () untuk dinamika lateral -x dan vertikal -y.

 lat = (0,2,3,5); alt = (1,4); #    Ax = (A[lat, :])[:, lat]; Bx = B[lat, 0]; Cx = C[0, lat]; Dx = D[0, 0]; Ay = (A[alt, :])[:, alt]; By = B[alt, 1]; Cy = C[1, alt]; Dy = D[1, 1]; #      K Qx1 = diag([1, 1, 1, 1, 1, 1]); Qu1a = diag([1, 1]); X,K,E =lqr(A, B, Qx1, Qu1a) K1a = matrix(K); #      x H1ax = ss(Ax - Bx*K1a[0,lat], Bx*K1a[0,lat]*xd[lat,:], Cx, Dx); (Yx, Tx) = step(H1ax, T=linspace(0,10,100)); #     y H1ay = ss(Ay - By*K1a[1,alt], By*K1a[1,alt]*yd[alt,:], Cy, Dy); (Yy, Ty) = step(H1ay, T=linspace(0,10,100)); 

Grafik ditampilkan secara berurutan, satu untuk setiap tugas.

 figure(); title("   x  y") plot(Tx.T, Yx.T, '-', Ty.T, Yy.T, '--'); plot([0, 10], [1, 1], 'k-'); axis([0, 10, -0.1, 1.4]); ylabel(' '); xlabel(''); legend(('x', 'y'), loc='lower right'); grid() show() 

Bagan:




Bagan:




Bagan:




Bagan:

Kesimpulan:

Implementasi optimasi LQR dalam pesawat yang disajikan dalam publikasi, dengan mempertimbangkan penggunaan fungsi lqr () baru, akan memungkinkan kita untuk meninggalkan penggunaan program berlisensi MATLAB ketika mempelajari bagian yang sesuai dari teori kontrol.

Referensi

1. Matematika pada jari: pengontrol linear-kuadratik.

2. Ruang keadaan dalam masalah merancang sistem kontrol yang optimal.

3. Menggunakan Pustaka Sistem Kontrol Python untuk merancang sistem kontrol otomatis.

4. Perpustakaan Sistem Kontrol Python.

5. Python - tidak dapat mengimpor modul slycot ke spyder (RuntimeError & ImportError).

6. Perpustakaan Sistem Kontrol Python.

7. Kriteria untuk kontrol optimal dan optimalisasi lqr dalam drive listrik.