Lyhyt dokumentaatio ohjelmasta, joka ratkaisee numeerisesti kahden pisteen reuna-arvotehtävän.
Tämä dokumentti jupyter-notebookina
Ongelma, joka tuottaa tässä ratkaistavan kahden-pisteen reuna-arvotehtävän
Ratkaisen reuna-arvotehtävän scikits.bvp_solver
-algoritmilla. Minun koneeseeni se installoitui paikkaan, josta python ei sitä löytänyt, joten tein allaolevan tempun. (Installoin komennolla pip3.6 install scikits.bvp_solver)
# -*- coding: utf-8 -*-
import sys
import os
# sys.path.append('/home/puudeli/0Pilvi/webPy/pylib/bvp_solver')
webpy = os.environ.get('WEBPY')
sys.path.append(webpy + 'pylib/bvp_solver')
# import pylib.bvp_solver
from solver import solve
from ProblemDefinition import ProblemDefinition
from Solution import Solution
from template_generator import get_template
import numpy
import time
from pylab import ioff
# from math import sin, cos
from timeit import default_timer as timer
from opt_dfuns import odefun, bcfun, init_guess # , gradodefun
from opt_params import par
from OptPiirto import opt_nosturi_piirto, tallenna_ohjaus
Numeerisen algoritmin konvergointi eli se, löytääkö algoritmi ratkaisun, riippuu alkuarvauksesta ja tehtävään liittyvistä parametreista. 'Käsityönä' koodaamani alkuarvaus antoi ratkaisun vain sellaisilla optimikustannuksen painokertoimilla, jotka eivät mielestäni olleet mielenkiintoisia. Niinpä tein funktion guess
, joka muuttaa optimoinnin tuloksen uudeksi alkuarvaukseksi bvp_solverin
toivomaan muotoon seuraavaa optimointia varten. Uusi alkuarvaus 'kestää' mielenkiintoisempia painokertoimia.
# guess palauttaa tilavektorin arvon ajanhetkellä t
def guess(t):
(tmax, dim) = guess_arr.shape
it = int((tmax - 1.0) / par.T_f * t)
return guess_arr[it, :]
Alempana määritelty funktio ratkaisu
hakee optimiohjauksen kutsumalla funktiota scikits.bvp_solver.solve
def solve(bvp_problem,
solution_guess,
initial_mesh = None,
parameter_guess = None,
max_subintervals = 300,
singular_term = None,
tolerance = 1.0e-6,
method = 4,
trace = 0,
error_on_fail = True):
"""Attempts to solve the supplied boundary value problem starting from the user supplied guess for the solution using BVP_SOLVER.
Parameters
----------
bvp_problem : :class:`ProblemDefinition`
Defines the boundary value problem to be solved.
solution_guess : :class:`Solution`, constant, array of values or function
A guess for the solution.
initial_mesh : castable to floating point ndarray
Points on the x-axis to use for the supplied solution guess, default is 10 evenly spaced points. Must not be supplied if solution_guess is a :class:`Solution` object.
parameter_guess : castable to floating point ndarray, shape (num_parameters)
Guesses for the unknown parameters. Must not be supplied if solution_guess is a :class:`Solution` object.
max_subintervals : int
Maximum number of points on the mesh before an error is returned.
singular_term : castable to floating point ndarray, shape(num_ODE, num_ODE)
Matrix that defines the singular term for the problem if one exist.
tolerance : positive float
Tolerance for size of defect of approximate solution.
method : {2, 4, 6}
Order of Runge-Kutta to use.
trace : {0, 1, 2}
Indicates verbosity of output. 0 for no output, 1 for some output, 2 for full output.
error_on_fail : logical
Indicates whether an exception should be raised if solving fails.
Returns
-------
sol : :class:`Solution`
Approximate problem solution.
Raises
------
ValueError
If bvp_problem failed validation in some way.
ValueError
If solving fails.
"""
def ratkaisu(i, s_paiva):
global guess_arr
start = timer()
solution = solve(
problem,
solution_guess=guess,
tolerance=1.0e-5,
method=6,
trace=2,
max_subintervals=500)
# x: aika-akseli, jolla on plotdim pistettä alku- ja loppuhetken välillä
x = numpy.linspace(problem.boundary_points[0], problem.boundary_points[1],
plotdim)
# y: ratkaisun mukaiset arvot tilamuuttujille kussakin aikapisteessä
y = solution(x)
# ratkaisusta uusi alkuarvaus
guess_arr = y.transpose()
end = timer()
# piirretään muuttujien kuvaajat ajan funktiona
opt_nosturi_piirto(s_paiva, i, x, y, plotdim)
# välitulostus, ratkaisun iterointiin kulunut aika ja
# muutaman mielenkiintoisen parametrin arvo
s1 = '\niteration: {0:2d} time: {1:.4g} '.format(i, end - start)
s2 = ' C_W {0:.4g}; C_F1 {1:.4g}; C_F2 {2:.4g}'
s2 = s2.format(par.C_W, par.C_F1, par.C_F2)
print(s1 + s2)
return x, y
Pythonilla kuvaajia voi piirrellä interaktiivisessa moodissa, jolloin grafiikkaohjelmisto jää odottamaan käyttäjältä lisäkomentoja. Minusta ohjelma "panttasi" joidenkin kuvaajien piirtämistä, joten käänsin interaktiivisen moodin pois päältä
# Interaktiivinen grafiikka pois käytöstä.
ioff()
plotdim
on aikapisteiden lukumäärä, joissa ratkaisu esitetään. Tämä on turhan paljon, mutta videopeli päivittää tilansa 48 kertaa sekunnissa, joten näin saan kätevästi automaattiohjauksen ohjausvoimat suoraan kullekin aikahetkellä.
plotdim = int(48.0 * par.T_f)
'Käsinkoodaamani' varsinainen alkuarvaus
guess_arr = numpy.array(
[init_guess(float(i) * par.T_f / 100.0) for i in range(0, 100)])
# kahden pisteen reuna-arvotehtävän määrittely
problem = ProblemDefinition(
num_ODE=4 * par.sysdim,
num_parameters=0,
num_left_boundary_conditions=2 * par.sysdim,
boundary_points=(0, par.T_f),
function=odefun,
# function_derivative=gradodefun,
boundary_conditions=bcfun)
class ProblemDefinition:
"""Defines a boundary value problem.
"""
def __init__(self,
num_ODE,
num_parameters,
num_left_boundary_conditions,
boundary_points,
function,
boundary_conditions,
function_derivative = None,
boundary_conditions_derivative = None):
"""
Parameters
----------
num_ODE : int
Number of first order ordinary differential equations in the problem.
num_parameters : int
Number of unknown parameters in the problem.
num_left_boundary_conditions : int
Number of boundary conditions enforced on the left boundary.
boundary_points : arraylike, shape(2)
Array that defines the two boundary points on the x axis.
function : function (see definition below)
A function which calculates the value of the ODE equations.
function (X, Y[, P]):
Parameters:
X : float
scalar value of x at which to evaluate the ODEs
Y : ndarray, shape(num_ODE)
current value of all variables
P : ndarray, shape(num_parameters)
value of all unknown parameters (only included if num_parameters > 0)
Returns:
ODE : ndarray, shape(num_ODE)
array of all ODEs
boundary_conditions : function (see definition below)
A function which calculates the difference between the boundary conditions and the actual variables currently calculated.
boundary_conditions(YA, YB[, P]):
Parameters:
YA : ndarray, shape(num_ODE)
value of all variables at the left boundary
YB : ndarray, shape(num_ODE)
value of all variables at the right boundary
P : ndarray, shape(num_parameters)
value of all unknown parameters (only used if num_parameters > 0)
Returns:
BCA : ndarray, shape(num_left_boundary_conditions)
difference between the boundary condition and variables at the left boundary
BCB : ndarray, shape(num_ODE + num_parameters - num_left_boundary_conditions)
array of the difference between the boundary condition and variables at the right boundary
function_derivative : optional function (see definition below)
A function which returns the partial derivatives of the function argument.
function (X, Y[, P]):
Parameters:
X : float
scalar value of x at which to evaluate the ODEs
Y : ndarray, shape(num_ODE)
current value of all variables
P : ndarray, shape(num_parameters)
value of all unknown parameters (only included if num_parameters > 0)
Returns:
dODE : ndarray, shape(num_ODE, num_ODE)
array of partial derivative of all ODEs with respect to all variables; index of ODEs is first, index of variables is second
dOdP : ndarray, shape(num_ODE, num_parameters)
array of partial derivative of all ODEs with respect to all unknown parameters; index of ODEs is first, index of parameters is second must not be returned if the problem does not include unknown parameters
boundary_conditions_drivative : optional function (see definition below)
A function which returns the partial derivatives of the boundary_conditions argument.
boundary_conditions(YA, YB[, P]):
Parameters:
YA : ndarray, shape(num_ODE)
value of all variables at the left boundary
YB : ndarray, shape(num_ODE)
value of all variables at the right boundary
P : ndarray, shape(num_parameters)
value of all unknown parameters (only used if num_parameters > 0)
Returns:
dBCA : ndarray, shape(num_left_boundary_conditions, num_ODE)
partial derivatives of the difference between the left boundary condition and the actual variables at the left boundary; boundary condition index is first and variable index is second
dBCB : ndarray, shape(num_ODE + num_parameters - num_left_boundary_conditions, num_ODE)
partial derivatives of the difference between the right boundary condition and the actual variables at the right boundary; boundary condition index is first and variable index is second
dBPA : ndarray, shape(num_left_boundary_conditions, num_parameters)
partial derivatives of the difference between the left boundary condition and the unknown parameters; boundary condition index is first and parameter index is second
dBPB : ndarray, shape(num_ODE + num_parameters - num_left_boundary_conditions, num_parameters)
partial derivatives of the difference between the right boundary condition and the unknown parameters; boundary condition index is first and parameter index is second
"""
paiva = time.localtime(time.time())
s_paiva = str(paiva.tm_year) + str(paiva.tm_mon) + str(paiva.tm_mday) + str(
paiva.tm_hour) + str(paiva.tm_min)
maxiter = 50
last = 0
# Muutetaan optimointikriteerin parametreja ja lasketaan uusi ratkaisu
# käyttäen edellistä ratkaisua alkuarvauksena
# Alla lisätään iteraatio iteraatiolta energian käytön painotusta.
# Voisi olla mielenkiintoista muutella ilmanvastusta ja kitkoja.
try:
for i in range(0, maxiter):
par.C_F1 = 0.8 * par.C_F1
par.C_F2 = 0.8 * par.C_F2
x, y = ratkaisu(i, s_paiva)
last = i
tallenna_ohjaus(x, y, plotdim)
finally:
print("se siitä")
# try: .. finally: rakenne
# Jos optimiratkaisua ei löydy, `bvp_solver` ilmoittaa suoritusvirheestä,
# jonka seurauksena suoritus siirtyy kohtaan finally:
iteration: 0 time: 10.25 C_W 1; C_F1 8; C_F2 1.6
iteration: 1 time: 10.03 C_W 1; C_F1 6.4; C_F2 1.28
iteration: 2 time: 6.639 C_W 1; C_F1 5.12; C_F2 1.024
iteration: 3 time: 9.691 C_W 1; C_F1 4.096; C_F2 0.8192
iteration: 4 time: 9.912 C_W 1; C_F1 3.277; C_F2 0.6554
iteration: 5 time: 7.844 C_W 1; C_F1 2.621; C_F2 0.5243
iteration: 6 time: 7.98 C_W 1; C_F1 2.097; C_F2 0.4194
iteration: 7 time: 6.208 C_W 1; C_F1 1.678; C_F2 0.3355
iteration: 8 time: 6.131 C_W 1; C_F1 1.342; C_F2 0.2684
iteration: 9 time: 7.279 C_W 1; C_F1 1.074; C_F2 0.2147
iteration: 10 time: 4.688 C_W 1; C_F1 0.859; C_F2 0.1718
iteration: 11 time: 5.815 C_W 1; C_F1 0.6872; C_F2 0.1374
iteration: 12 time: 4.436 C_W 1; C_F1 0.5498; C_F2 0.11
iteration: 13 time: 4.354 C_W 1; C_F1 0.4398; C_F2 0.08796
iteration: 14 time: 5.581 C_W 1; C_F1 0.3518; C_F2 0.07037
iteration: 15 time: 4.376 C_W 1; C_F1 0.2815; C_F2 0.05629
iteration: 16 time: 4.391 C_W 1; C_F1 0.2252; C_F2 0.04504
iteration: 17 time: 4.393 C_W 1; C_F1 0.1801; C_F2 0.03603
iteration: 18 time: 4.441 C_W 1; C_F1 0.1441; C_F2 0.02882
iteration: 19 time: 6.045 C_W 1; C_F1 0.1153; C_F2 0.02306
iteration: 20 time: 5.901 C_W 1; C_F1 0.09223; C_F2 0.01845
iteration: 21 time: 6.255 C_W 1; C_F1 0.07379; C_F2 0.01476
iteration: 22 time: 6.553 C_W 1; C_F1 0.05903; C_F2 0.01181
iteration: 23 time: 6.79 C_W 1; C_F1 0.04722; C_F2 0.009445
iteration: 24 time: 7.638 C_W 1; C_F1 0.03778; C_F2 0.007556
iteration: 25 time: 8.318 C_W 1; C_F1 0.03022; C_F2 0.006045
iteration: 26 time: 7.508 C_W 1; C_F1 0.02418; C_F2 0.004836
iteration: 27 time: 133.4 C_W 1; C_F1 0.01934; C_F2 0.003869 se siitä
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-10-e2b7f9898b60> in <module> 13 par.C_F1 = 0.8 * par.C_F1 14 par.C_F2 = 0.8 * par.C_F2 ---> 15 x, y = ratkaisu(i, s_paiva) 16 last = i 17 tallenna_ohjaus(x, y, plotdim) <ipython-input-5-2857bc3d3672> in ratkaisu(i, s_paiva) 2 global guess_arr 3 start = timer() ----> 4 solution = solve( 5 problem, 6 solution_guess=guess, ~/0Pilvi/webPy/pylib/bvp_solver/solver.py in solve(bvp_problem, solution_guess, initial_mesh, parameter_guess, max_subintervals, singular_term, tolerance, method, trace, error_on_fail) 172 # check to see if there was a problem with the solution 173 if error_on_fail and calculatedSolution._successIndicator == -1: --> 174 raise ValueError("Boundary value problem solving failed. Run with trace = 1 or 2 for more information.") 175 176 return calculatedSolution ValueError: Boundary value problem solving failed. Run with trace = 1 or 2 for more information.