SciPy (ออกเสียงว่า sai pie) เป็นแพ็คเกจคณิตศาสตร์แบบ numpy ซึ่งรวมถึงไลบรารี C และ Fortran ด้วย SciPy เปลี่ยนเซสชัน Python แบบโต้ตอบของคุณให้เป็นสภาพแวดล้อมวิทยาศาสตร์ข้อมูลที่สมบูรณ์ เช่น MATLAB, IDL, Octave, R หรือ SciLab
ในบทความนี้ เราจะดูเทคนิคพื้นฐานของการเขียนโปรแกรมทางคณิตศาสตร์ - การแก้ปัญหาการปรับให้เหมาะสมตามเงื่อนไขสำหรับฟังก์ชันสเกลาร์ของตัวแปรหลายตัวโดยใช้แพ็คเกจ scipy.optimize อัลกอริธึมการเพิ่มประสิทธิภาพที่ไม่มีข้อจำกัดได้ถูกกล่าวถึงแล้ว . สามารถรับความช่วยเหลือโดยละเอียดและอัปเดตเพิ่มเติมเกี่ยวกับฟังก์ชัน scipy ได้เสมอโดยใช้คำสั่ง help(), Shift+Tab หรือใน .
การแนะนำ
อินเทอร์เฟซทั่วไปสำหรับการแก้ปัญหาการปรับให้เหมาะสมทั้งแบบมีเงื่อนไขและแบบไม่จำกัดในแพ็คเกจ scipy.optimize มีให้โดยฟังก์ชัน minimize(). อย่างไรก็ตามเป็นที่ทราบกันดีว่าไม่มีวิธีการสากลในการแก้ปัญหาทั้งหมดดังนั้นการเลือกวิธีการที่เหมาะสมจึงตกอยู่บนไหล่ของผู้วิจัยเช่นเคย
อัลกอริทึมการปรับให้เหมาะสมที่เหมาะสมถูกระบุโดยใช้อาร์กิวเมนต์ของฟังก์ชัน minimize(..., method="").
สำหรับการเพิ่มประสิทธิภาพตามเงื่อนไขของฟังก์ชันของตัวแปรหลายตัว มีวิธีการใช้งานดังต่อไปนี้:
trust-constr— ค้นหาค่าต่ำสุดในท้องถิ่นในภูมิภาคความเชื่อมั่น , ;SLSQP— การโปรแกรมกำลังสองตามลำดับที่มีข้อจำกัด วิธีนิวตันสำหรับการแก้ระบบลากรองจ์ .TNC- มีข้อจำกัดของนิวตันที่ถูกตัดทอน จำนวนการวนซ้ำที่จำกัด เหมาะสำหรับฟังก์ชันไม่เชิงเส้นที่มีตัวแปรอิสระจำนวนมาก .L-BFGS-B- วิธีการจากทีม Broyden–Fletcher–Goldfarb–Shanno ซึ่งนำไปใช้โดยลดการใช้หน่วยความจำเนื่องจากการโหลดเวกเตอร์บางส่วนจากเมทริกซ์ Hessian , .COBYLA— การเพิ่มประสิทธิภาพที่จำกัดของ MARE โดยการประมาณเชิงเส้น การเพิ่มประสิทธิภาพที่จำกัดด้วยการประมาณเชิงเส้น (โดยไม่มีการคำนวณการไล่ระดับสี) .
ขึ้นอยู่กับวิธีการที่เลือก เงื่อนไขและข้อจำกัดในการแก้ปัญหาจะถูกตั้งค่าแตกต่างกันไป:
- วัตถุคลาส
Boundsสำหรับวิธีการ L-BFGS-B, TNC, SLSQP, trust-constr; - รายการ
(min, max)สำหรับวิธีการเดียวกัน L-BFGS-B, TNC, SLSQP, trust-constr; - วัตถุหรือรายการวัตถุ
LinearConstraint,NonlinearConstraintสำหรับ COBYLA, SLSQP, วิธี trust-constr; - พจนานุกรมหรือรายการพจนานุกรม
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}สำหรับวิธี COBYLA, SLSQP
โครงร่างบทความ:
1) พิจารณาการใช้อัลกอริธึมการปรับให้เหมาะสมตามเงื่อนไขในภูมิภาคที่เชื่อถือได้ (method=”trust-constr”) โดยมีข้อจำกัดที่ระบุเป็นอ็อบเจ็กต์ Bounds, LinearConstraint, NonlinearConstraint ;
2) พิจารณาการเขียนโปรแกรมตามลำดับโดยใช้วิธีกำลังสองน้อยที่สุด (method = "SLSQP") โดยมีข้อจำกัดที่ระบุไว้ในรูปแบบของพจนานุกรม {'type', 'fun', 'jac', 'args'};
3) วิเคราะห์ตัวอย่างการปรับให้เหมาะสมของผลิตภัณฑ์ที่ผลิตโดยใช้ตัวอย่างของเว็บสตูดิโอ
วิธีการปรับให้เหมาะสมตามเงื่อนไข = "trust-constr"
การนำวิธีการไปใช้ trust-constr ขึ้นอยู่กับ สำหรับปัญหาเกี่ยวกับข้อจำกัดของรูปแบบความเสมอภาคและต่อ สำหรับปัญหาข้อจำกัดในรูปของความไม่เท่าเทียมกัน ทั้งสองวิธีถูกนำมาใช้โดยอัลกอริธึมเพื่อค้นหาค่าต่ำสุดในพื้นที่ในภูมิภาคความเชื่อมั่น และเหมาะสมอย่างยิ่งสำหรับปัญหาขนาดใหญ่
การกำหนดทางคณิตศาสตร์ของปัญหาการหาค่าต่ำสุดในรูปแบบทั่วไป:
สำหรับข้อจำกัดความเท่าเทียมกันที่เข้มงวด ขอบเขตล่างจะถูกตั้งค่าให้เท่ากับขอบเขตบน .
สำหรับข้อจำกัดทางเดียว ขีดจำกัดบนหรือล่างจะถูกตั้งค่าไว้ np.inf โดยมีเครื่องหมายตรงกัน
ปล่อยให้จำเป็นต้องค้นหาค่าต่ำสุดของฟังก์ชัน Rosenbrock ที่รู้จักของตัวแปรสองตัว:
ในกรณีนี้ มีการตั้งค่าข้อจำกัดต่อไปนี้ในโดเมนของคำจำกัดความ:
ในกรณีของเรา มีวิธีแก้ปัญหาเฉพาะจุดตรงประเด็น ซึ่งใช้ได้เฉพาะข้อจำกัดข้อแรกและข้อที่สี่เท่านั้น
มาดูข้อจำกัดจากล่างขึ้นบนแล้วดูว่าเราจะเขียนมันในรูปแบบ scipy ได้อย่างไร
ข้อ จำกัด и มากำหนดมันโดยใช้วัตถุ Bounds
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])ข้อ จำกัด и ลองเขียนมันในรูปแบบเชิงเส้น:
มากำหนดข้อจำกัดเหล่านี้เป็นวัตถุ LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])และสุดท้ายข้อจำกัดไม่เชิงเส้นในรูปแบบเมทริกซ์:
เรากำหนดเมทริกซ์จาโคเบียนสำหรับข้อจำกัดนี้และผลรวมเชิงเส้นของเมทริกซ์ Hessian ด้วยเวกเตอร์ใดๆ :
ตอนนี้เราสามารถกำหนดข้อจำกัดไม่เชิงเส้นเป็นวัตถุได้ 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)หากมีขนาดใหญ่ คุณสามารถระบุเมทริกซ์ในรูปแบบกระจัดกระจายได้:
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)หรือเป็นวัตถุ 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)เมื่อคำนวณเมทริกซ์ Hessian ต้องใช้ความพยายามมาก คุณสามารถใช้คลาสได้ . มีกลยุทธ์ต่อไปนี้: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hessian สามารถคำนวณได้โดยใช้ผลต่างอันจำกัด:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')เมทริกซ์จาโคเบียนสำหรับข้อจำกัดสามารถคำนวณได้โดยใช้ผลต่างอันจำกัด อย่างไรก็ตาม ในกรณีนี้ ไม่สามารถคำนวณเมทริกซ์ Hessian โดยใช้ผลต่างอันจำกัดได้ Hessian ต้องถูกกำหนดให้เป็นฟังก์ชันหรือใช้คลาส HessianUpdateStrategy
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())วิธีแก้ไขปัญหาการปรับให้เหมาะสมมีลักษณะดังนี้:
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]หากจำเป็น สามารถกำหนดฟังก์ชันสำหรับการคำนวณ Hessian ได้โดยใช้คลาส 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)หรือผลคูณของ Hessian และเวกเตอร์ใดๆ ผ่านพารามิเตอร์ 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)อีกวิธีหนึ่ง สามารถประมาณอนุพันธ์ตัวแรกและตัวที่สองของฟังก์ชันที่กำลังปรับให้เหมาะสมได้ ตัวอย่างเช่น Hessian สามารถประมาณได้โดยใช้ฟังก์ชัน SR1 (การประมาณกึ่งนิวตัน) การไล่ระดับสีสามารถประมาณได้ด้วยผลต่างอันจำกัด
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)วิธีการปรับให้เหมาะสมตามเงื่อนไข = "SLSQP"
วิธี SLSQP ได้รับการออกแบบมาเพื่อแก้ไขปัญหาการลดฟังก์ชันให้เหลือน้อยที่สุดในรูปแบบ:
ที่ไหน и — ชุดดัชนีของสำนวนที่อธิบายข้อจำกัดในรูปแบบของความเสมอภาคหรือความไม่เท่าเทียมกัน — ชุดของขอบเขตล่างและบนสำหรับโดเมนของคำจำกัดความของฟังก์ชัน
ข้อจำกัดเชิงเส้นและไม่เชิงเส้นอธิบายไว้ในรูปแบบของพจนานุกรมพร้อมปุ่มต่างๆ type, fun и 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])
}การค้นหาขั้นต่ำจะดำเนินการดังนี้:
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 ]ตัวอย่างการเพิ่มประสิทธิภาพ
ในการเชื่อมต่อกับการเปลี่ยนไปใช้โครงสร้างทางเทคโนโลยีที่ห้า เรามาดูการเพิ่มประสิทธิภาพการผลิตโดยใช้ตัวอย่างของสตูดิโอเว็บ ซึ่งทำให้เรามีรายได้เพียงเล็กน้อยแต่มั่นคง ลองจินตนาการว่าเราเป็นผู้อำนวยการห้องครัวที่ผลิตผลิตภัณฑ์สามประเภท:
- x0 - ขายหน้า Landing Page จาก 10 tr.
- x1 - เว็บไซต์องค์กร ตั้งแต่ 20 tr.
- x2 - ร้านค้าออนไลน์จาก 30 tr.
ทีมงานที่เป็นมิตรของเราประกอบด้วยรุ่นน้องสี่คน คนกลางสองคน และผู้อาวุโสหนึ่งคน กองทุนเวลาทำงานรายเดือน:
- มิถุนายน:
4 * 150 = 600 чел * час, - ตรงกลาง:
2 * 150 = 300 чел * час, - อาวุโส:
150 чел * час.
ปล่อยให้ผู้เยาว์คนแรกใช้เวลา (0, 1, 2) ชั่วโมงในการพัฒนาและปรับใช้ไซต์ประเภทหนึ่ง (x10, x20, x30), กลาง - (7, 15, 20), อาวุโส - (5, 10, 15 ) ชั่วโมงแห่งช่วงเวลาที่ดีที่สุดในชีวิตของคุณ
เช่นเดียวกับกรรมการทั่วไป เราต้องการเพิ่มผลกำไรรายเดือนให้สูงสุด ขั้นตอนแรกสู่ความสำเร็จคือการเขียนฟังก์ชันวัตถุประสงค์ value เป็นจำนวนรายได้จากผลิตภัณฑ์ที่ผลิตต่อเดือน:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]นี่ไม่ใช่ข้อผิดพลาด เมื่อค้นหาค่าสูงสุด ฟังก์ชันวัตถุประสงค์จะถูกย่อให้เล็กสุดด้วยเครื่องหมายตรงกันข้าม
ขั้นตอนต่อไปคือการห้ามไม่ให้พนักงานของเราทำงานหนักเกินไปและจำกัดเวลาทำงาน:
สิ่งที่เทียบเท่า:
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]])
}ข้อจำกัดอย่างเป็นทางการคือผลผลิตของผลิตภัณฑ์ต้องเป็นค่าบวกเท่านั้น:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])และสุดท้าย ข้อสันนิษฐานที่เป็นสีดอกกุหลาบที่สุดคือ เนื่องจากราคาที่ต่ำและมีคุณภาพสูง ลูกค้าที่พึงพอใจจึงต่อคิวรอเราอยู่ตลอดเวลา เราสามารถเลือกปริมาณการผลิตรายเดือนได้ด้วยตนเอง โดยอิงจากการแก้ปัญหาการเพิ่มประสิทธิภาพที่มีข้อจำกัด 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]ปัดเศษให้เป็นจำนวนเต็มแล้วคำนวณปริมาณพายในแต่ละเดือนด้วยการกระจายผลิตภัณฑ์ที่เหมาะสมที่สุด x = (8, 6, 3) :
- มิถุนายน:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - ตรงกลาง:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - อาวุโส:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
สรุป: เพื่อให้ผู้อำนวยการได้รับสูงสุดที่สมควรจะเป็นวิธีที่ดีที่สุดในการสร้างแลนดิ้งเพจ 8 หน้า ไซต์ขนาดกลาง 6 แห่ง และร้านค้า 3 แห่งต่อเดือน ในกรณีนี้ผู้อาวุโสจะต้องไถโดยไม่เงยหน้าขึ้นจากเครื่อง โหลดตรงกลางจะอยู่ที่ประมาณ 2/3 ส่วนรุ่นน้องน้อยกว่าครึ่งหนึ่ง
ข้อสรุป
บทความนี้สรุปเทคนิคพื้นฐานในการทำงานกับแพ็คเกจ scipy.optimizeใช้เพื่อแก้ปัญหาการย่อเล็กสุดแบบมีเงื่อนไข ส่วนตัวผมใช้ scipy มีจุดประสงค์เพื่อการศึกษาเท่านั้น ซึ่งเป็นเหตุให้ตัวอย่างที่ให้มามีลักษณะเป็นการ์ตูน
มีทฤษฎีและตัวอย่างเสมือนจริงมากมายในหนังสือของ I.L. Akulich เรื่อง "การเขียนโปรแกรมทางคณิตศาสตร์ในตัวอย่างและปัญหา" แอปพลิเคชั่นที่ฮาร์ดคอร์มากขึ้น scipy.optimize เพื่อสร้างโครงสร้าง 3 มิติจากชุดรูปภาพ () สามารถดูได้ใน .
แหล่งข้อมูลหลักคือ ผู้ที่ต้องการมีส่วนร่วมในการแปลในส่วนนี้และส่วนอื่นๆ scipy ยินดีต้อนรับสู่ .
ขอบคุณ เพื่อมีส่วนร่วมในการจัดทำสิ่งพิมพ์
ที่มา: will.com