Commit 9dbb0241 authored by Jorn Baayen's avatar Jorn Baayen

transcribe() working again

parent 3837bc50
from casadi import MX, Function, jacobian, vertcat, reshape, mtimes, substitute, interpolant
from casadi import MX, Function, jacobian, vertcat, reshape, mtimes, substitute, interpolant, transpose, repmat
import numpy as np
import logging
......@@ -10,7 +10,7 @@ def is_affine(e, v):
return (f.sparsity_jac(0, 0).nnz() == 0)
def nullvertcat(L):
def nullvertcat(*L):
"""
Like vertcat, but creates an MX with consistent dimensions even if L is empty.
"""
......@@ -26,8 +26,8 @@ def reduce_matvec(e, v):
This reduces the number of nodes required to represent the linear operations.
"""
Af = Function("Af", [], [jacobian(e, v)])
A = Af()[0]
Af = Function("Af", [MX()], [jacobian(e, v)])
A = Af(MX())
return reshape(mtimes(A, v), e.shape)
......@@ -43,8 +43,10 @@ def reduce_matvec_plus_b(e, v):
def interpolate(ts, xs, t, equidistant, mode=0):
if mode == 0:
return interpolant(ts, xs, t, equidistant)
if False: # TODO mode == 0:
print(ts)
print(xs)
return interpolant('interpolant', 'linear', [ts], xs, {'lookup_mode': 'exact'})(t)
else:
if mode == 1:
xs = xs[:-1] # block-forward
......@@ -54,8 +56,8 @@ def interpolate(ts, xs, t, equidistant, mode=0):
if t.size1() > 1:
t_ = MX.sym('t')
xs_ = MX.sym('xs', xs.size1())
f = Function('interpolant', [t_, xs_], [mul(transpose((t_ >= ts[:-1]) * (t_ < ts[1:])), xs_)])
f = f.map('interpolant_map', t.size1())
return transpose(f([transpose(t), repmat(xs, 1, t.size1())])[0])
f = Function('interpolant', [t_, xs_], [mtimes(transpose((t_ >= ts[:-1]) * (t_ < ts[1:])), xs_)])
f = f.map(t.size1(), 'serial')
return transpose(f(transpose(t), repmat(xs, 1, t.size1())))
else:
return mul(transpose((t >= ts[:-1]) * (t < ts[1:])), xs)
\ No newline at end of file
return mtimes(transpose((t >= ts[:-1]) * (t < ts[1:])), xs)
\ No newline at end of file
......@@ -60,6 +60,9 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
for var in itertools.chain(self.dae_variables['states'], self.dae_variables['algebraics'], self.dae_variables['control_inputs'], self.dae_variables['constant_inputs'], self.dae_variables['parameters'], self.dae_variables['time']):
self._variables[var.name()] = var
# Call super
super().__init__(**kwargs)
@abstractmethod
def times(self, variable=None):
"""
......@@ -354,14 +357,14 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
# Aggregate ensemble data
ensemble_aggregate = {}
ensemble_aggregate["parameters"] = horzcat(nullvertcat(*l) for l in ensemble_parameter_values)
ensemble_aggregate["initial_constant_inputs"] = horzcat(nullvertcat(*[float(d["constant_inputs"][variable.name()][0])
for variable in self.dae_variables['constant_inputs']]) for d in ensemble_store)
ensemble_aggregate["initial_state"] = horzcat(d["initial_state"] for d in ensemble_store)
ensemble_aggregate["parameters"] = horzcat(*[nullvertcat(*l) for l in ensemble_parameter_values])
ensemble_aggregate["initial_constant_inputs"] = horzcat(*[nullvertcat(*[float(d["constant_inputs"][variable.name()][0])
for variable in self.dae_variables['constant_inputs']]) for d in ensemble_store])
ensemble_aggregate["initial_state"] = horzcat(*[d["initial_state"] for d in ensemble_store])
ensemble_aggregate["initial_state"] = reduce_matvec(ensemble_aggregate["initial_state"], self.solver_input)
ensemble_aggregate["initial_derivatives"] = horzcat(d["initial_derivatives"] for d in ensemble_store)
ensemble_aggregate["initial_derivatives"] = horzcat(*[d["initial_derivatives"] for d in ensemble_store])
ensemble_aggregate["initial_derivatives"] = reduce_matvec(ensemble_aggregate["initial_derivatives"], self.solver_input)
ensemble_aggregate["initial_path_variables"] = horzcat(d["initial_path_variables"] for d in ensemble_store)
ensemble_aggregate["initial_path_variables"] = horzcat(*[d["initial_path_variables"] for d in ensemble_store])
ensemble_aggregate["initial_path_variables"] = reduce_matvec(ensemble_aggregate["initial_path_variables"], self.solver_input)
if (self._dae_residual_function_collocated is None) and (self._integrator_step_function is None):
......@@ -411,7 +414,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
dae_residual_collocated = vertcat(*dae_residual_collocated)
# Check linearity of collocated part
if self.check_collocation_linearity and dae_residual_collocated.size1() > 0:
if False: #TODO self.check_collocation_linearity and dae_residual_collocated.size1() > 0:
# Check linearity of collocation constraints, which is a necessary condition for the optimization problem to be convex
self.linear_collocation = True
if not is_affine(dae_residual_collocated, vertcat(*
......@@ -446,7 +449,6 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
dae_residual_function_integrated = Function('dae_residual_function_integrated', [I, I0, ensemble_parameters, vertcat(*[C0[i] for i in range(len(collocated_variables))] + [CI0[i] for i in range(len(
self.dae_variables['constant_inputs']))] + [dt_sym] + collocated_variables + collocated_derivatives + self.dae_variables['constant_inputs'] + self.dae_variables['time'])], [dae_residual_integrated], function_options)
dae_residual_function_integrated = dae_residual_function_integrated.expand()
options = self.integrator_options()
self._integrator_step_function = rootfinder('integrator_step_function', 'newton', dae_residual_function_integrated, options)
......@@ -455,13 +457,12 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
if len(collocated_variables) > 0:
self._dae_residual_function_collocated = Function('dae_residual_function_collocated', [ensemble_parameters, vertcat(*
integrated_variables + collocated_variables + integrated_derivatives + collocated_derivatives + self.dae_variables['constant_inputs'] + self.dae_variables['time'])], [dae_residual_collocated], function_options)
self._dae_residual_function_collocated = self._dae_residual_function_collocated.expand()
if len(integrated_variables) > 0:
integrator_step_function = self._integrator_step_function
if len(collocated_variables) > 0:
dae_residual_function_collocated = self._dae_residual_function_collocated
dae_residual_collocated_size = dae_residual_function_collocated.outputExpr(0).size1()
dae_residual_collocated_size = dae_residual_function_collocated.mx_out(0).size1()
else:
dae_residual_collocated_size = 0
......@@ -555,7 +556,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
if len(collocated_variables) > 0:
if theta < 1:
# Obtain state vector
[dae_residual_0] = dae_residual_function_collocated([ensemble_parameters,
[dae_residual_0] = dae_residual_function_collocated.call([ensemble_parameters,
vertcat(integrated_states_0,
collocated_states_0,
integrated_finite_differences,
......@@ -565,7 +566,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
False, True)
if theta > 0:
# Obtain state vector
[dae_residual_1] = dae_residual_function_collocated([ensemble_parameters,
[dae_residual_1] = dae_residual_function_collocated.call([ensemble_parameters,
vertcat(integrated_states_1,
collocated_states_1,
integrated_finite_differences,
......@@ -581,7 +582,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
accumulated_Y.append(
(1 - theta) * dae_residual_0 + theta * dae_residual_1)
accumulated_Y.extend(path_objective_function([ensemble_parameters,
accumulated_Y.extend(path_objective_function.call([ensemble_parameters,
vertcat(integrated_states_1,
collocated_states_1,
integrated_finite_differences,
......@@ -591,7 +592,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
path_variables_1)],
False, True))
accumulated_Y.extend(path_constraints_function([ensemble_parameters,
accumulated_Y.extend(path_constraints_function.call([ensemble_parameters,
vertcat(integrated_states_1,
collocated_states_1,
integrated_finite_differences,
......@@ -608,13 +609,11 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
[accumulated_X, accumulated_U, ensemble_parameters], [vertcat(*accumulated_Y)], function_options)
if len(integrated_variables) > 0:
accumulation = accumulated.mapaccum(
'accumulation', n_collocation_times - 1)
accumulation = accumulated.mapsum(n_collocation_times - 1)
else:
# Fully collocated problem. Use map(), so that we can use
# parallelization along the time axis.
accumulation = accumulated.map(
'accumulation', n_collocation_times - 1, {'parallelization': 'openmp'})
accumulation = accumulated.map(n_collocation_times - 1, 'openmp')
# Start collecting constraints
f = []
......@@ -627,10 +626,9 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
initial_residual_with_params_fun = Function('initial_residual', [ensemble_parameters, vertcat(*self.dae_variables['states'] + self.dae_variables['algebraics'] + self.dae_variables[
'control_inputs'] + integrated_derivatives + collocated_derivatives + self.dae_variables['constant_inputs'] + self.dae_variables['time'])], [vertcat(*[dae_residual, initial_residual])],
function_options)
initial_residual_with_params_fun = initial_residual_with_params_fun.expand()
self._initial_residual_with_params_fun_map = initial_residual_with_params_fun.map('initial_residual_with_params_fun_map', self.ensemble_size)
self._initial_residual_with_params_fun_map = initial_residual_with_params_fun.map(self.ensemble_size)
initial_residual_with_params_fun_map = self._initial_residual_with_params_fun_map
[res] = initial_residual_with_params_fun_map([ensemble_aggregate["parameters"], vertcat(*[ensemble_aggregate["initial_state"], ensemble_aggregate["initial_derivatives"], ensemble_aggregate["initial_constant_inputs"], repmat([0.0], 1, self.ensemble_size)])], False, True)
[res] = initial_residual_with_params_fun_map.call([ensemble_aggregate["parameters"], vertcat(*[ensemble_aggregate["initial_state"], ensemble_aggregate["initial_derivatives"], ensemble_aggregate["initial_constant_inputs"], repmat([0.0], 1, self.ensemble_size)])], False, True)
res = vec(res)
g.append(res)
......@@ -746,8 +744,8 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
# Map to all time steps
logger.info("Mapping")
[integrators_and_collocation_and_path_constraints] = accumulation(
[accumulation_X0, accumulation_U, repmat(parameters, 1, n_collocation_times - 1)])
integrators_and_collocation_and_path_constraints = accumulation(
accumulation_X0, accumulation_U, repmat(parameters, 1, n_collocation_times - 1))
if integrators_and_collocation_and_path_constraints.size2() > 0:
integrators = integrators_and_collocation_and_path_constraints[:len(integrated_variables), :]
collocation_constraints = vec(integrators_and_collocation_and_path_constraints[len(integrated_variables):len(
......@@ -776,8 +774,6 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
# Add collocation constraints
if collocation_constraints.size1() > 0:
#if self._affine_collocation_constraints:
# collocation_constraints = reduce_matvec_plus_b(collocation_constraints, self.solver_input)
g.append(collocation_constraints)
zeros = np.zeros(collocation_constraints.size1())
lbg.extend(zeros)
......@@ -835,7 +831,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
# Objective
f_member = self.objective(ensemble_member)
if path_objective.size1() > 0:
initial_path_objective = path_objective_function([parameters,
initial_path_objective = path_objective_function.call([parameters,
vertcat(initial_state
, initial_derivatives
, initial_constant_inputs,
......@@ -869,7 +865,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
path_constraints = self.path_constraints(ensemble_member)
if len(path_constraints) > 0:
# We need to evaluate the path constraints at t0, as the initial time is not included in the accumulation.
[initial_path_constraints] = path_constraints_function([parameters,
[initial_path_constraints] = path_constraints_function.call([parameters,
vertcat(initial_state
, initial_derivatives
, initial_constant_inputs,
......@@ -1274,7 +1270,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
# Discretization parameters
control_size = self._control_size
ensemble_member_size = self._state_size / self.ensemble_size
ensemble_member_size = int(self._state_size / self.ensemble_size)
# Extract control inputs
results = {}
......@@ -1354,7 +1350,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
# Look up transcribe_problem() state.
X = self.solver_input
control_size = self._control_size
ensemble_member_size = self._state_size / self.ensemble_size
ensemble_member_size = int(self._state_size / self.ensemble_size)
# Find state vector
vector = None
......@@ -1395,7 +1391,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
t0 = self.initial_time
X = self.solver_input
control_size = self._control_size
ensemble_member_size = self._state_size / self.ensemble_size
ensemble_member_size = int(self._state_size / self.ensemble_size)
# Fetch appropriate symbol, or value.
canonical, sign = self.alias_relation.canonical_signed(variable)
......@@ -1496,7 +1492,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
# Look up transcribe_problem() state.
X = self.solver_input
control_size = self._control_size
ensemble_member_size = self._state_size / self.ensemble_size
ensemble_member_size = int(self._state_size / self.ensemble_size)
# Compute position in state vector
offset = control_size + ensemble_member * ensemble_member_size
......@@ -1637,7 +1633,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
# We have a special symbol for t0 derivatives
X = self.solver_input
control_size = self._control_size
ensemble_member_size = self._state_size / self.ensemble_size
ensemble_member_size = int(self._state_size / self.ensemble_size)
canonical, sign = self.alias_relation.canonical_signed(variable)
try:
......@@ -1690,7 +1686,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass = A
f = Function('f', [vertcat(*states_and_path_variables), vertcat(*derivatives),
vertcat(*self.dae_variables['constant_inputs']), vertcat(*self.dae_variables['parameters']),
self.dae_variables['time'][0]], [expr])
fmap = f.map('fmap', len(self.times()))
fmap = f.map(len(self.times()))
# Discretization settings
collocation_times = self.times()
......
# cython: embedsignature=True
from casadi import MX, substitute, repmat, vertcat, depends_on
from casadi import MX, substitute, repmat, vertcat, depends_on, veccat
from pymola.backends.casadi.api import transfer_model
from collections import OrderedDict
import numpy as np
......@@ -72,7 +72,7 @@ class ModelicaMixin(OptimizationProblem):
else:
if v.symbol.name() in kwargs.get('lookup_tables', []):
self._mx['lookup_tables'].append(v.symbol)
elif v.fixed == False:
elif v.fixed == True:
self._mx['constant_inputs'].append(v.symbol)
else:
self._mx['control_inputs'].append(v.symbol)
......@@ -106,12 +106,13 @@ class ModelicaMixin(OptimizationProblem):
v.symbol.name(), alias))
# Initialize nominals and types
# TODO needs initialized metadata
self._nominals = {}
self._discrete = {}
for v in itertools.chain(self._pymola_model.states, self._pymola_model.alg_states, self._pymola_model.inputs):
sym_name = v.symbol.name()
if v.nominal != 0 and v.nominal != 1:
self._nominals[sym_name] = abs(float(nominal))
self._nominals[sym_name] = abs(float(v.nominal))
if logger.getEffectiveLevel() == logging.DEBUG:
logger.debug("ModelicaMixin: Set nominal value for variable {} to {}".format(
......@@ -119,6 +120,19 @@ class ModelicaMixin(OptimizationProblem):
self._discrete[sym_name] = v.python_type != float
# Initialize dae and initial residuals
inputs = [v.symbol for v in self._pymola_model.inputs]
self._dae_residual = self._pymola_model.dae_residual_function(self._mx['time'][0],
veccat(*self._mx['states']), veccat(*self._mx['derivatives']), veccat(*self._mx['algebraics']), veccat(*inputs), MX(), veccat(*self._mx['parameters']))
if self._dae_residual is None:
self._dae_residual = MX()
self._initial_residual = self._pymola_model.initial_residual_function(self._mx['time'][0],
veccat(*self._mx['states']), veccat(*self._mx['derivatives']), veccat(*self._mx['algebraics']), veccat(*inputs), MX(), veccat(*self._mx['parameters']))
if self._initial_residual is None:
self._initial_residual = MX()
# Call parent class first for default behaviour.
super(ModelicaMixin, self).__init__(**kwargs)
......@@ -138,10 +152,10 @@ class ModelicaMixin(OptimizationProblem):
# Eliminate constant symbols from model, replacing them with the values
# specified in the model.
compiler_options['replace_constants'] = True
compiler_options['replace_constant_values'] = True
# Replace any parameter expressions into the model.
compiler_options['replace_parameter_expressions'] = 'True'
compiler_options['replace_parameter_expressions'] = True
# Eliminate variables starting with underscores.
compiler_options['eliminable_variable_expression'] = r'_\w+'
......@@ -150,8 +164,7 @@ class ModelicaMixin(OptimizationProblem):
compiler_options['detect_aliases'] = True
# Cache the model on disk
# TODO set to False to runt into unset attributes
compiler_options['cache'] = False
compiler_options['cache'] = True
# Done
return compiler_options
......@@ -211,7 +224,7 @@ class ModelicaMixin(OptimizationProblem):
# Initial conditions obtained from start attributes.
for v in self._pymola_model.states:
# TODO nan values
initial_state[sym.name()] = sym.start
initial_state[v.symbol.name()] = v.start
return initial_state
......@@ -238,7 +251,7 @@ class ModelicaMixin(OptimizationProblem):
if M is None:
M = 1
m_ = v.min
m_ = MX(v.min)
if not m_.is_constant():
[m] = substitute([m_], self._mx['parameters'], parameter_values)
if m.is_constant():
......@@ -248,7 +261,7 @@ class ModelicaMixin(OptimizationProblem):
if np.isfinite(m_): # TODO vector values
m = m_
M_ = v.max
M_ = MX(v.max)
if not M_.is_constant():
[M] = substitute([M_], self._mx['parameters'], parameter_values)
if M.is_constant():
......
......@@ -72,6 +72,9 @@ class OptimizationProblem(metaclass = ABCMeta):
Base class for all optimization problems.
"""
def __init__(self, **kwargs):
self._mixed_integer = False
def optimize(self, preprocessing=True, postprocessing=True, log_solver_failure_as_error=True):
"""
Perform one initialize-transcribe-solve-finalize cycle.
......@@ -111,14 +114,10 @@ class OptimizationProblem(metaclass = ABCMeta):
my_solver = options['solver']
del options['solver']
# Expand option
expand = options['expand']
del options['expand']
# Already consumed
del options['optimized_num_dir']
nlpsol_options = {'expand': expand, my_solver: options}
nlpsol_options = {my_solver: options}
if self._mixed_integer:
nlpsol_options['discrete'] = discrete
......
......@@ -11,7 +11,7 @@ model TestModel
output Real z;
input Real x_delayed(fixed=true);
input Real x_delayed(fixed=false);
output Real switched;
......
......@@ -12,7 +12,7 @@ model TestModelWithInitial
output Real z;
input Real x_delayed(fixed=true);
input Real x_delayed(fixed=false);
output Real switched;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment