SciPy, optimisasi dengan kondisi

SciPy (diucapkan sai pie) adalah paket matematika berbasis numpy yang juga mencakup perpustakaan C dan Fortran. Dengan SciPy, sesi Python interaktif berubah menjadi lingkungan pemrosesan data yang lengkap seperti MATLAB, IDL, Octave, R, atau SciLab.

Pada artikel ini, kami akan mempertimbangkan teknik dasar pemrograman matematika - memecahkan masalah optimisasi bersyarat untuk fungsi skalar dari beberapa variabel menggunakan paket scipy.optimize. Algoritma optimasi tanpa syarat sudah dipertimbangkan dalam artikel terakhir . Bantuan yang lebih terperinci dan terkini tentang fungsi-fungsi yang membingungkan selalu dapat diperoleh dengan menggunakan perintah help (), Shift + Tab, atau dalam dokumentasi resmi .

Pendahuluan

Antarmuka umum untuk memecahkan masalah optimasi bersyarat dan tanpa syarat dalam paket scipy.optimize disediakan oleh fungsi minimize() . Namun, diketahui bahwa metode universal untuk menyelesaikan semua masalah tidak ada, oleh karena itu, seperti biasa, pilihan metode yang memadai ada pada peneliti.

Algoritme pengoptimalan yang sesuai ditentukan menggunakan argumen ke fungsi minimize(..., method="") .

Untuk optimisasi bersyarat dari fungsi beberapa variabel, metode berikut tersedia:

• trust-constr - mencari minimum lokal di bidang kepercayaan. Artikel wiki ; Artikel Habr ;
• SLSQP - pemrograman kuadratik berurutan dengan kendala, metode Newton untuk memecahkan sistem Lagrange. Artikel wiki .
• TNC - Truncated Newton Constrained, sejumlah iterasi terbatas, bagus untuk fungsi nonlinier dengan sejumlah besar variabel independen. Artikel wiki .
• L-BFGS-B - metode dari empat Broyden - Fletcher - Goldfarb - Shanno, diimplementasikan dengan konsumsi memori yang berkurang karena pemuatan sebagian vektor dari matriks Hessian. Artikel wiki , artikel Habr .
• COBYLA - MARES Optimalisasi Terkini Dengan Pendekatan Linear, optimasi terbatas dengan pendekatan linier (tanpa perhitungan gradien). Artikel wiki .

Bergantung pada metode yang dipilih, kondisi dan batasan untuk menyelesaikan masalah diatur secara berbeda:

• objek kelas Bounds untuk metode L-BFGS-B, TNC, SLSQP, trust-constr;
• daftar (min, max) untuk metode yang sama L-BFGS-B, TNC, SLSQP, trust-constr;
• objek atau daftar objek LinearConstraint , NonlinearConstraint untuk metode COBYLA, SLSQP, trust-constr;
• kamus atau daftar kamus {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} untuk metode COBYLA, SLSQP.

Garis besar artikel:
1) Pertimbangkan penerapan algoritma optimasi bersyarat di area kepercayaan (metode = "kepercayaan-konstr") dengan batasan yang ditentukan dalam bentuk Bounds , LinearConstraint , objek NonlinearConstraint ;
2) Pertimbangkan pemrograman kuadrat terkecil berurutan (metode = "SLSQP") dengan batasan yang ditentukan dalam kamus {'type', 'fun', 'jac', 'args'} ;
3) Untuk membongkar contoh optimasi produk pada contoh studio web.

Metode optimasi bersyarat = "trust-constr"

Implementasi metode trust-constr didasarkan pada EQSQP untuk masalah dengan pembatasan bentuk kesetaraan dan pada TRIP untuk masalah dengan pembatasan dalam bentuk ketidaksetaraan. Kedua metode diimplementasikan oleh algoritma untuk menemukan minimum lokal di bidang kepercayaan dan cocok untuk tugas skala besar.

Formulasi matematis dari masalah pencarian minimum dalam bentuk umum:

$$

$$\ min_x f (x)$$

$$

$$c ^ l \ leq c (x) \ leq c ^ u,$$

$$

$$x ^ l \ leq x \ leq x ^ u$$

Untuk kendala kesetaraan yang ketat, batas bawah ditetapkan sama dengan batas atas $$$c ^ l_j = c ^ u_j$ .
Untuk pembatasan satu arah, batas atas atau bawah ditetapkan oleh np.inf dengan tanda yang sesuai.
Biarkan diperlukan untuk menemukan fungsi Rosenbrock minimum dari dua variabel yang diketahui:

$$

$$\ min_ {x_0, x_1} 100 (x_0 - x_1 ^ 2) ^ 2 + (1-x_0) ^ 2$$

Dalam hal ini, batasan berikut ditetapkan pada domain definisinya:

$$

$$x_0 ^ 2 + x_1 \ leq 1$$

$$

$$x_0 ^ 2 - x_1 \ leq 1$$

$$

$$2x_0 + x_1 = 1$$

$$

$$x_0 + 2x_1 \ leq 1$$

$$

$$0 \ leq x_0 \ leq 1$$

$$

$$-0.5 \ leq x_1 \ leq $ 2.$$

Dalam kasus kami, ada solusi unik pada saat itu $$$[x_0, x_1] = [0,4149, 0,1701]$ hanya pembatasan pertama dan keempat yang valid.
Mari kita pergi melalui pembatasan dari bawah ke atas dan mempertimbangkan bagaimana mereka dapat ditulis dalam scipy.
Keterbatasan $$$0 \ leq x_0 \ leq 1$ dan $$$-0.5 \ leq x_1 \ leq $ 2.$ didefinisikan menggunakan objek Batas.

 from scipy.optimize import Bounds bounds = Bounds ([0, -0.5], [1.0, 2.0]) 

Keterbatasan $$$x_0 + 2 x_1 \ leq 1$ dan $$$2 x_0 + x_1 = 1$ kami menulis dalam bentuk linear:

$$

$$\ begin {bmatrix} - \ infty \\ 1 \ end {bmatrix} \ leq \ begin {bmatrix} 1 & 2 \\ 2 & 1 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \ end {bmatrix } \ leq \ begin {bmatrix} 1 \\ 1 \ end {bmatrix}$$

Tentukan batasan ini sebagai objek LinearConstraint:

 import numpy as np from scipy.optimize import LinearConstraint linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1]) 

Dan akhirnya, batasan nonlinear dalam bentuk matriks:

$$

$$c (x) = \ begin {bmatrix} x_0 ^ 2 + x_1 \\ x_0 ^ 2 - x_1 \ end {bmatrix} \ leq \ begin {bmatrix} 1 \\ 1 \ end {bmatrix}$$

Kami mendefinisikan matriks Jacobi untuk pembatasan ini dan kombinasi linear dari matriks Hessian dengan vektor sewenang-wenang $$$v$ :

$$

$$J (x) = \ begin {bmatrix} 2x_0 & 1 \\ 2x_0 & -1 \ end {bmatrix}$$

$$

$$H (x, v) = \ jumlah \ limit_ {i = 0} ^ 1 v_i \ nabla ^ 2 c_i (x) = v_0 \ begin {bmatrix} 2 & 0 \\ 2 & 0 \ end {bmatrix} + v_1 \ begin {bmatrix} 2 & 0 \\ 2 & 0 \ end {bmatrix}$$

Sekarang kita dapat mendefinisikan batasan nonlinear sebagai objek NonlinearConstraint :

 from scipy.optimize import NonlinearConstraint def cons_f(x): return [x[0]**2 + x[1], x[0]**2 - x[1]] def cons_J(x): return [[2*x[0], 1], [2*x[0], -1]] def cons_H(x, v): return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]]) nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H) 

Jika ukurannya besar, matriks juga dapat ditentukan dalam bentuk yang jarang:

 from scipy.sparse import csc_matrix def cons_H_sparse(x, v): return v[0]*csc_matrix([[2, 0], [0, 0]]) + v[1]*csc_matrix([[2, 0], [0, 0]]) nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H_sparse) 

atau sebagai objek LinearOperator :

 from scipy.sparse.linalg import LinearOperator def cons_H_linear_operator(x, v): def matvec(p): return np.array([p[0]*2*(v[0]+v[1]), 0]) return LinearOperator((2, 2), matvec=matvec) nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H_linear_operator) 

Saat perhitungan matriks Hessian $$$H (x, v)$ mahal, Anda dapat menggunakan kelas HessianUpdateStrategy . Strategi berikut tersedia: BFGS dan SR1 .

 from scipy.optimize import BFGS nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS()) 

Goni juga dapat dihitung menggunakan perbedaan hingga:

 nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point') 

Matriks Jacobi untuk batasan juga dapat dihitung menggunakan perbedaan hingga. Namun, dalam hal ini, matriks Hessian tidak dapat dihitung menggunakan perbedaan hingga. Hessian harus didefinisikan sebagai fungsi atau menggunakan kelas HessianUpdateStrategy.

 nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ()) 

Solusi untuk masalah optimasi adalah sebagai berikut:

 from scipy.optimize import minimize from scipy.optimize import rosen, rosen_der, rosen_hess, rosen_hess_prod x0 = np.array([0.5, 0]) res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess, constraints=[linear_constraint, nonlinear_constraint], options={'verbose': 1}, bounds=bounds) print(res.x) 

 `gtol` termination condition is satisfied. Number of iterations: 12, function evaluations: 8, CG iterations: 7, optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 0.033 s. [0.41494531 0.17010937] 

Jika perlu, fungsi perhitungan Goni dapat didefinisikan menggunakan kelas LinearOperator

 def rosen_hess_linop(x): def matvec(p): return rosen_hess_prod(x, p) return LinearOperator((2, 2), matvec=matvec) res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess_linop, constraints=[linear_constraint, nonlinear_constraint], options={'verbose': 1}, bounds=bounds) print(res.x) 

atau produk dari Hessian dan vektor sewenang-wenang melalui parameter hessp :

 res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hessp=rosen_hess_prod, constraints=[linear_constraint, nonlinear_constraint], options={'verbose': 1}, bounds=bounds) print(res.x) 

Atau, turunan pertama dan kedua dari fungsi yang dioptimalkan dapat dihitung kira-kira. Misalnya, Hessian dapat diperkirakan menggunakan fungsi SR1 (pendekatan kuasi-Newtonian). Gradien dapat diperkirakan dengan perbedaan hingga.

 from scipy.optimize import SR1 res = minimize(rosen, x0, method='trust-constr', jac="2-point", hess=SR1(), constraints=[linear_constraint, nonlinear_constraint], options={'verbose': 1}, bounds=bounds) print(res.x) 

Metode optimasi bersyarat = "SLSQP"

Metode SLSQP dirancang untuk memecahkan masalah meminimalkan fungsi dalam formulir

$$

$$\ min_x f (x)$$

$$

$$c_j (x) = 0, j \ in \ mathcal {E}$$

$$

$$c_j (x) \ geq 0, j \ in \ mathcal {I}$$

$$

$$lb_i \ leq x_i \ leq ub_i, i = 1, ..., N$$

Dimana $$$\ mathcal {E}$ dan $$$\ mathcal {I}$ - set indeks ekspresi yang menggambarkan kendala dalam bentuk persamaan atau ketidaksetaraan $$$[lb_i, ub_i]$ - set batas bawah dan atas untuk domain definisi fungsi.

Batasan linear dan nonlinier dijelaskan dalam bentuk kamus dengan type kunci, fun dan jac .

 ineq_cons = {'type': 'ineq', 'fun': lambda x: np.array ([1 - x [0] - 2 * x [1], 1 - x [0] ** 2 - x [1], 1 - x [0] ** 2 + x [1]]), 'jac': lambda x: np.array ([[- 1.0, -2.0], [-2 * x [0], -1.0], [-2 * x [0], 1.0]]) } eq_cons = {'type': 'eq', 'fun': lambda x: np.array ([2 * x [0] + x [1] - 1]), 'jac': lambda x: np.array ([2.0, 1.0]) } 

Pencarian minimum dilakukan sebagai berikut:

 x0 = np.array([0.5, 0]) res = minimize(rosen, x0, method='SLSQP', jac=rosen_der, constraints=[eq_cons, ineq_cons], options={'ftol': 1e-9, 'disp': True}, bounds=bounds) print(res.x) 

 Optimization terminated successfully. (Exit mode 0) Current function value: 0.34271757499419825 Iterations: 4 Function evaluations: 5 Gradient evaluations: 4 [0.41494475 0.1701105 ] 

Contoh optimasi

Sehubungan dengan transisi ke struktur teknologi kelima, kami akan mempertimbangkan optimalisasi produksi menggunakan contoh web studio, yang memberi kami penghasilan kecil namun stabil. Mari kita perkenalkan diri kita sebagai direktur galeri, yang menghasilkan tiga jenis produk:

• x0 - jual pendaratan, mulai 10 tr
• x1 - situs perusahaan, mulai 20 tr
• x2 - toko online, mulai 30 tr

Tim kerja kami yang ramah meliputi empat Juns, dua Middle dan satu Senior. Dana waktu kerja mereka selama sebulan:

• Juni: 4 * 150 = 600 * ,
• middles: 2 * 150 = 300 * ,
• Senior: 150 *

Misalkan junior pertama yang dihabiskan untuk mengembangkan dan menggunakan satu situs tipe (x0, x1, x2) harus menghabiskan (10, 20, 30) jam, menengah - (7, 15, 20), senior - (5, 10, 15) jam terbaik waktu hidupnya.

Seperti direktur normal lainnya, kami ingin memaksimalkan laba bulanan kami. Langkah pertama menuju sukses adalah menuliskan value fungsi tujuan sebagai jumlah pendapatan dari produksi bulanan:

 def value(x): return - 10*x[0] - 20*x[1] - 30*x[2] 

Ini bukan kesalahan, ketika mencari maksimum, fungsi objektif diminimalkan dengan tanda sebaliknya.

Langkah selanjutnya adalah melarang pemrosesan untuk karyawan kami dan memperkenalkan batasan pada dana waktu kerja:

$$

$$\ begin {bmatrix} 10 & 20 & 30 \\ 7 & 15 & 20 \\ 5 & 10 & 15 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \\ x_2 \ end {bmatrix} \ leq \ mulai {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix}$$

Yang setara:

$$

$$\ begin {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix} - \ begin {bmatrix} 10 & 20 & 30 \\ 7 & 15 & 20 \\ 5 & 10 & 15 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \\ x_2 \ end {bmatrix} \ geq 0$$

 ineq_cons = {'type': 'ineq', 'fun': lambda x: np.array ([600 - 10 * x [0] - 20 * x [1] - 30 * x[2], 300 - 7 * x [0] - 15 * x [1] - 20 * x[2], 150 - 5 * x [0] - 10 * x [1] - 15 * x[2]]) } 

Pembatasan formal - hasilnya hanya positif:

 bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf]) 

Dan akhirnya, asumsi paling optimis - karena harga rendah dan kualitas tinggi, sederetan pelanggan yang puas terus berbaris untuk kami. Kami dapat memilih sendiri volume produksi bulanan, berdasarkan solusi masalah optimisasi bersyarat dengan scipy.optimize :

 x0 = np.array([10, 10, 10]) res = minimize(value, x0, method='SLSQP', constraints=ineq_cons, bounds=bnds) print(res.x) 

 [7.85714286 5.71428571 3.57142857] 

Kami mengitari lollipop ke bilangan bulat terdekat dan menghitung beban pendayung bulanan pada distribusi produksi yang optimal x = (8, 6, 3) :

• Juni: 8 * 10 + 6 * 20 + 3 * 30 = 290 * ;
• middles: 8 * 7 + 6 * 15 + 3 * 20 = 206 * ;
• Senior: 8 * 5 + 6 * 10 + 3 * 15 = 145 * .

Kesimpulan: bagi direktur untuk mendapatkan nilai maksimum yang layak, optimal untuk membuat 8 pendaratan, 6 situs menengah dan 3 toko per bulan. Pada saat yang sama, seigneur harus membajak tanpa melepaskan diri dari mesin, beban tengah sekitar 2/3, kurang dari setengah Juni.

Kesimpulan

Artikel ini menguraikan teknik dasar untuk bekerja dengan paket scipy.optimize yang digunakan untuk memecahkan masalah minimalisasi bersyarat. Secara pribadi, saya menggunakan scipy murni untuk tujuan akademis, jadi contoh ini sangat lucu.

Banyak teori dan contoh winrar dapat ditemukan, misalnya, dalam buku I. L. Akulich "Pemrograman Matematika dalam Contoh dan Masalah". Aplikasi yang lebih hardcore dari scipy.optimize untuk membangun struktur 3D untuk satu set gambar ( artikel di hub ) dapat ditemukan di buku masak-scipy .

Sumber utama informasi adalah docs.scipy.org. Jika Anda ingin menerjemahkan ini dan bagian lain dari scipy ke dalam terjemahan, selamat datang di GitHub .

Terima kasih kepada mephistopheies yang berkontribusi pada publikasi ini.