# ___________________________________________________________________________
#
# Pyomo: Python Optimization Modeling Objects
# Copyright (c) 2008-2025
# National Technology and Engineering Solutions of Sandia, LLC
# Under the terms of Contract DE-NA0003525 with National Technology and
# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
# rights in this software.
# This software is distributed under the 3-clause BSD License.
# ___________________________________________________________________________
import logging
import math
# If the user has numpy then the collocation points and the a matrix for
# the Runge-Kutta basis formulation will be calculated as needed.
# If the user does not have numpy then these values will be read from a
# stored dictionary for up to 10 collocation points.
from pyomo.common.dependencies import numpy, numpy_available
from pyomo.common.collections import ComponentSet
from pyomo.core.base import Transformation, TransformationFactory
from pyomo.core import Var, ConstraintList, Expression, Objective
from pyomo.dae import ContinuousSet, DerivativeVar, Integral
from pyomo.dae.misc import generate_finite_elements
from pyomo.dae.misc import generate_colloc_points
from pyomo.dae.misc import expand_components
from pyomo.dae.misc import create_partial_expression
from pyomo.dae.misc import add_discretization_equations
from pyomo.dae.misc import add_continuity_equations
from pyomo.dae.misc import block_fully_discretized
from pyomo.dae.misc import get_index_information
from pyomo.dae.diffvar import DAE_Error
from pyomo.common.config import ConfigBlock, ConfigValue, PositiveInt, In
logger = logging.getLogger('pyomo.dae')
def _lagrange_radau_transform(v, s):
ncp = s.get_discretization_info()['ncp']
adot = s.get_discretization_info()['adot']
def _fun(i):
tmp = list(s)
idx = s.ord(i) - 1
if idx == 0: # Don't apply this equation at initial point
raise IndexError("list index out of range")
low = s.get_lower_element_boundary(i)
lowidx = s.ord(low) - 1
return sum(
v(tmp[lowidx + j])
* adot[j][idx - lowidx]
* (1.0 / (tmp[lowidx + ncp] - tmp[lowidx]))
for j in range(ncp + 1)
)
return _fun
def _lagrange_radau_transform_order2(v, s):
ncp = s.get_discretization_info()['ncp']
adotdot = s.get_discretization_info()['adotdot']
def _fun(i):
tmp = list(s)
idx = s.ord(i) - 1
if idx == 0:
# Don't apply this equation at initial point
raise IndexError("list index out of range")
low = s.get_lower_element_boundary(i)
lowidx = s.ord(low) - 1
return sum(
v(tmp[lowidx + j])
* adotdot[j][idx - lowidx]
* (1.0 / (tmp[lowidx + ncp] - tmp[lowidx]) ** 2)
for j in range(ncp + 1)
)
return _fun
def _lagrange_legendre_transform(v, s):
ncp = s.get_discretization_info()['ncp']
adot = s.get_discretization_info()['adot']
def _fun(i):
tmp = list(s)
idx = s.ord(i) - 1
if idx == 0:
# Don't apply this equation at initial point
raise IndexError("list index out of range")
elif i in s.get_finite_elements():
# Don't apply at finite element points continuity equations
# added later
raise IndexError("list index out of range")
low = s.get_lower_element_boundary(i)
lowidx = s.ord(low) - 1
return sum(
v(tmp[lowidx + j])
* adot[j][idx - lowidx]
* (1.0 / (tmp[lowidx + ncp + 1] - tmp[lowidx]))
for j in range(ncp + 1)
)
return _fun
def _lagrange_legendre_transform_order2(v, s):
ncp = s.get_discretization_info()['ncp']
adotdot = s.get_discretization_info()['adotdot']
def _fun(i):
tmp = list(s)
idx = s.ord(i) - 1
if idx == 0:
# Don't apply this equation at initial point
raise IndexError("list index out of range")
elif i in s.get_finite_elements():
# Don't apply at finite element points continuity equations
# added later
raise IndexError("list index out of range")
low = s.get_lower_element_boundary(i)
lowidx = s.ord(low) - 1
return sum(
v(tmp[lowidx + j])
* adotdot[j][idx - lowidx]
* (1.0 / (tmp[lowidx + ncp + 1] - tmp[lowidx]) ** 2)
for j in range(ncp + 1)
)
return _fun
[docs]
def conv(a, b):
if len(a) == 0 or len(b) == 0:
raise ValueError("Cannot convolve an empty list")
ans = []
m = len(a)
n = len(b)
for k in range(m + n - 1):
val = 0
j = max(0, k - n)
stop = min(k, m)
while j <= stop:
if j < m and (k - j) < n:
val += a[j] * b[k - j]
j += 1
ans.insert(k, val)
return ans
[docs]
def calc_cp(alpha, beta, k):
gamma = []
factorial = math.factorial
for i in range(k + 1):
num = factorial(alpha + k) * factorial(alpha + beta + k + i)
denom = factorial(alpha + i) * factorial(k - i) * factorial(i)
gamma.insert(i, num / denom)
poly = []
for i in range(k + 1):
if i == 0:
poly.insert(i, gamma[i])
else:
prod = [1]
j = 1
while j <= i:
prod = conv(prod, [1, -1])
j += 1
while len(poly) < len(prod):
poly.insert(0, 0)
prod = [gamma[i] * t for t in prod]
poly = [sum(pair) for pair in zip(poly, prod)]
cp = numpy.roots(poly)
return numpy.sort(cp).tolist()
# BLN: This is a legacy function that was used to calculate the collocation
# constants for an alternative form of the collocation equations described
# in Biegler's nonlinear programming book. The difference being whether the
# state or the derivative is approximated using lagrange polynomials. With
# the addition of PDE support and chained discretizations in Pyomo.DAE 2.0
# this function is no longer used but kept here for future reference.
#
# def calc_omega(cp):
# a = []
# for i in range(len(cp)):
# ptmp = []
# tmp = 0
# for j in range(len(cp)):
# if j != i:
# row = []
# row.insert(0, 1 / (cp[i] - cp[j]))
# row.insert(1, -cp[j] / (cp[i] - cp[j]))
# ptmp.insert(tmp, row)
# tmp += 1
# p = [1]
# for j in range(len(cp) - 1):
# p = conv(p, ptmp[j])
# pint = numpy.polyint(p)
# arow = []
# for j in range(len(cp)):
# arow.append(numpy.polyval(pint, cp[j]))
# a.append(arow)
# return a
[docs]
def calc_adot(cp, order=1):
a = []
for i in range(len(cp)):
ptmp = []
tmp = 0
for j in range(len(cp)):
if j != i:
row = []
row.insert(0, 1 / (cp[i] - cp[j]))
row.insert(1, -cp[j] / (cp[i] - cp[j]))
ptmp.insert(tmp, row)
tmp += 1
p = [1]
for j in range(len(cp) - 1):
p = conv(p, ptmp[j])
pder = numpy.polyder(p, order)
arow = []
for j in range(len(cp)):
arow.append(float(numpy.polyval(pder, cp[j])))
a.append(arow)
return a
[docs]
def calc_afinal(cp):
afinal = []
for i in range(len(cp)):
ptmp = []
tmp = 0
for j in range(len(cp)):
if j != i:
row = []
row.insert(0, 1 / (cp[i] - cp[j]))
row.insert(1, -cp[j] / (cp[i] - cp[j]))
ptmp.insert(tmp, row)
tmp += 1
p = [1]
for j in range(len(cp) - 1):
p = conv(p, ptmp[j])
afinal.append(float(numpy.polyval(p, 1.0)))
return afinal