Skip to content

Commit

Permalink
Add DiscreteTimeApproximation for converting general continuous-time …
Browse files Browse the repository at this point in the history
…systems
  • Loading branch information
wei-chen-li committed Feb 19, 2025
1 parent 36a6cc9 commit 6d44217
Show file tree
Hide file tree
Showing 8 changed files with 565 additions and 22 deletions.
1 change: 1 addition & 0 deletions bindings/pydrake/systems/BUILD.bazel
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,7 @@ drake_pybind_library(
cc_srcs = ["primitives_py.cc"],
package_info = PACKAGE_INFO,
py_deps = [
":analysis_py",
":framework_py",
":module_py",
"//bindings/pydrake:trajectories_py",
Expand Down
36 changes: 23 additions & 13 deletions bindings/pydrake/systems/primitives_py.cc
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
#include "drake/systems/primitives/constant_vector_source.h"
#include "drake/systems/primitives/demultiplexer.h"
#include "drake/systems/primitives/discrete_derivative.h"
#include "drake/systems/primitives/discrete_time_approximation.h"
#include "drake/systems/primitives/discrete_time_delay.h"
#include "drake/systems/primitives/discrete_time_integrator.h"
#include "drake/systems/primitives/first_order_low_pass_filter.h"
Expand Down Expand Up @@ -52,6 +53,7 @@ PYBIND11_MODULE(primitives, m) {
m.doc() = "Bindings for the primitives portion of the Systems framework.";
constexpr auto& doc = pydrake_doc.drake.systems;

py::module::import("pydrake.systems.analysis");
py::module::import("pydrake.systems.framework");
py::module::import("pydrake.trajectories");

Expand Down Expand Up @@ -240,6 +242,18 @@ PYBIND11_MODULE(primitives, m) {
py::arg("C") = Eigen::MatrixXd(), py::arg("D") = Eigen::MatrixXd(),
py::arg("time_period") = 0.0, doc.LinearSystem.ctor.doc_5args);

m.def("DiscreteTimeApproximation",
overload_cast_explicit<std::unique_ptr<LinearSystem<T>>,
const LinearSystem<T>&, double>(&DiscreteTimeApproximation),
py::arg("system"), py::arg("time_period"),
doc.DiscreteTimeApproximation.doc_linearsystem);

m.def("DiscreteTimeApproximation",
overload_cast_explicit<std::unique_ptr<AffineSystem<T>>,
const AffineSystem<T>&, double>(&DiscreteTimeApproximation),
py::arg("system"), py::arg("time_period"),
doc.DiscreteTimeApproximation.doc_affinesystem);

DefineTemplateClassWithDefault<MatrixGain<T>, LinearSystem<T>>(
m, "MatrixGain", GetPyParam<T>(), doc.MatrixGain.doc)
.def(py::init<const Eigen::Ref<const MatrixXd>&>(), py::arg("D"),
Expand Down Expand Up @@ -617,7 +631,7 @@ PYBIND11_MODULE(primitives, m) {
type_visit(bind_common_scalar_types, CommonScalarPack{});

// N.B. Capturing `&doc` should not be required; workaround per #9600.
auto bind_non_symbolic_scalar_types = [m, &doc](auto dummy) {
auto bind_non_symbolic_scalar_types = [&m, &doc](auto dummy) {
using T = decltype(dummy);

DefineTemplateClassWithDefault<LinearTransformDensity<T>, LeafSystem<T>>(m,
Expand Down Expand Up @@ -741,6 +755,14 @@ PYBIND11_MODULE(primitives, m) {
&TimeVaryingAffineSystem<T>::configure_random_state,
py::arg("covariance"),
doc.TimeVaryingAffineSystem.configure_random_state.doc);

m.def("DiscreteTimeApproximation",
overload_cast_explicit<std::unique_ptr<System<T>>, const System<T>&,
double, double, const drake::systems::SimulatorConfig&>(
&DiscreteTimeApproximation),
py::arg("system"), py::arg("time_period"), py::arg("time_offset") = 0.0,
py::arg("integrator_config") = drake::systems::SimulatorConfig(),
doc.DiscreteTimeApproximation.doc_system);
};
type_visit(bind_non_symbolic_scalar_types, NonSymbolicScalarPack{});

Expand Down Expand Up @@ -801,18 +823,6 @@ PYBIND11_MODULE(primitives, m) {

m.def("IsDetectable", &IsDetectable, py::arg("sys"),
py::arg("threshold") = std::nullopt, doc.IsDetectable.doc);

m.def("DiscreteTimeApproximation",
overload_cast_explicit<std::unique_ptr<LinearSystem<double>>,
const LinearSystem<double>&, double>(&DiscreteTimeApproximation),
py::arg("system"), py::arg("time_period"),
doc.DiscreteTimeApproximation.doc_linearsystem);

m.def("DiscreteTimeApproximation",
overload_cast_explicit<std::unique_ptr<AffineSystem<double>>,
const AffineSystem<double>&, double>(&DiscreteTimeApproximation),
py::arg("system"), py::arg("time_period"),
doc.DiscreteTimeApproximation.doc_affinesystem);
} // NOLINT(readability/fn_size)

} // namespace pydrake
Expand Down
18 changes: 15 additions & 3 deletions bindings/pydrake/systems/test/primitives_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -924,7 +924,8 @@ def test_vector_log_sink(self, T):
all(logger.GetLog(logger_context) == logger.FindLog(context)
for logger, logger_context in loggers_and_contexts))

def test_discrete_time_approximation(self):
@numpy_compare.check_all_types
def test_discrete_time_approximation_linearsystem(self, T):
A = np.array([[0, 1], [0, 0]])
B = np.array([0, 1])
f0 = np.array([2, 1])
Expand All @@ -943,15 +944,15 @@ def test_discrete_time_approximation(self):
def assert_array_close(x, y): np.testing.assert_allclose(
np.squeeze(x), np.squeeze(y), atol=1e-10)

continuous_system = LinearSystem(A, B, C, D)
continuous_system = LinearSystem_[T](A, B, C, D)
discrete_system = DiscreteTimeApproximation(
system=continuous_system, time_period=h)
assert_array_close(discrete_system.A(), Ad)
assert_array_close(discrete_system.B(), Bd)
assert_array_close(discrete_system.C(), Cd)
assert_array_close(discrete_system.D(), Dd)

continuous_system = AffineSystem(A, B, f0, C, D, y0)
continuous_system = AffineSystem_[T](A, B, f0, C, D, y0)
discrete_system = DiscreteTimeApproximation(
system=continuous_system, time_period=h)
assert_array_close(discrete_system.A(), Ad)
Expand All @@ -960,3 +961,14 @@ def assert_array_close(x, y): np.testing.assert_allclose(
assert_array_close(discrete_system.C(), Cd)
assert_array_close(discrete_system.D(), Dd)
assert_array_close(discrete_system.y0(), y0d)

def test_discrete_time_approximation_system(self):
sys = DiscreteTimeApproximation(
system=Integrator_[float](1),
time_period=0.01, time_offset=0.0)
sys.ToScalarType[AutoDiffXd]()

sys = DiscreteTimeApproximation(
system=Integrator_[AutoDiffXd](1),
time_period=0.01, time_offset=0.0)
sys.ToScalarType[float]()
27 changes: 27 additions & 0 deletions systems/primitives/BUILD.bazel
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ drake_cc_package_library(
":constant_vector_source",
":demultiplexer",
":discrete_derivative",
":discrete_time_approximation",
":discrete_time_delay",
":discrete_time_integrator",
":first_order_low_pass_filter",
Expand Down Expand Up @@ -121,6 +122,17 @@ drake_cc_library(
],
)

drake_cc_library(
name = "discrete_time_approximation",
srcs = ["discrete_time_approximation.cc"],
hdrs = ["discrete_time_approximation.h"],
deps = [
"//systems/analysis:simulator",
"//systems/analysis:simulator_config_functions",
"//systems/framework",
],
)

drake_cc_library(
name = "discrete_time_delay",
srcs = ["discrete_time_delay.cc"],
Expand Down Expand Up @@ -490,6 +502,21 @@ drake_cc_googletest(
],
)

drake_cc_googletest(
name = "discrete_time_approximation_test",
deps = [
":constant_vector_source",
":discrete_time_approximation",
":integrator",
":zero_order_hold",
"//multibody/parsing",
"//multibody/plant",
"//systems/analysis:simulator",
"//systems/controllers:linear_quadratic_regulator",
"//systems/framework",
],
)

drake_cc_googletest(
name = "discrete_time_delay_test",
deps = [
Expand Down
182 changes: 182 additions & 0 deletions systems/primitives/discrete_time_approximation.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,182 @@
#include "drake/systems/primitives/discrete_time_approximation.h"

#include <memory>

#include "drake/systems/analysis/simulator.h"
#include "drake/systems/analysis/simulator_config_functions.h"
#include "drake/systems/framework/leaf_system.h"

namespace drake {
namespace systems {

namespace {

template <typename T>
class DiscreteTimeSystem final : public LeafSystem<T> {
public:
DRAKE_NO_COPY_NO_MOVE_NO_ASSIGN(DiscreteTimeSystem);

DiscreteTimeSystem(std::unique_ptr<System<T>> system, double time_period,
double time_offset,
const SimulatorConfig& integrator_config)
: LeafSystem<T>(SystemTypeTag<DiscreteTimeSystem>{}),
continuous_system_(std::move(system)),
continuous_context_(continuous_system_->CreateDefaultContext()),
time_period_(time_period),
time_offset_(time_offset),
integrator_config_(integrator_config),
simulator__(Simulator(*continuous_system_)) {
this->Initialize();
}

void SetRandomParameters(const Context<T>& context, Parameters<T>* parameters,
RandomGenerator* generator) const override {
continuous_system_->SetRandomParameters(
*continuous_context_, &(continuous_context_->get_mutable_parameters()),
generator);
parameters->SetFrom(continuous_context_->get_parameters());
}

void SetRandomState(const Context<T>& context, State<T>* state,
RandomGenerator* generator) const override {
continuous_system_->SetRandomState(
*continuous_context_, &(continuous_context_->get_mutable_state()),
generator);
state->get_mutable_discrete_state().get_mutable_vector().SetFromVector(
continuous_context_->get_continuous_state_vector().CopyToVector());
}

// Scalar-converting copy constructor. See @ref system_scalar_conversion.
template <typename U>
explicit DiscreteTimeSystem(const DiscreteTimeSystem<U>& other)
: DiscreteTimeSystem<T>(
other.continuous_system_->template ToScalarType<T>(),
other.time_period_, other.time_offset_, other.integrator_config_) {}

private:
template <typename U>
friend class DiscreteTimeSystem;

void Initialize() {
// Configure integrator.
ApplySimulatorConfig(integrator_config_, &simulator__);
integrator_ = &(simulator__.get_mutable_integrator());
integrator_->reset_context(continuous_context_.get());
integrator_->Initialize();

// Set name.
const std::string& name = continuous_system_->get_name();
this->set_name(name.empty() ? "discrete-time approximation"
: "discrete-time approximated " + name);

// Declare input ports.
for (int i = 0; i < continuous_system_->num_input_ports(); ++i) {
const InputPort<T>& port = continuous_system_->get_input_port(i);
this->DeclareInputPort(port.get_name(), port.get_data_type(), port.size(),
port.get_random_type());
}

// Declare output ports.
for (int i = 0; i < continuous_system_->num_output_ports(); ++i) {
const OutputPort<T>& port = continuous_system_->get_output_port(i);
if (port.get_data_type() == PortDataType::kVectorValued) {
this->DeclareVectorOutputPort(
port.get_name(), port.size(),
[this, &port](const Context<T>& discrete_context,
BasicVector<T>* out) {
this->CopyDiscreteContextToContinuousContext(discrete_context);
out->SetFromVector(port.Eval(*(this->continuous_context_)));
});
} else if (port.get_data_type() == PortDataType::kAbstractValued) {
this->DeclareAbstractOutputPort(
port.get_name(),
[&port] {
return port.Allocate();
},
[this, &port](const Context<T>& discrete_context,
AbstractValue* out) {
this->CopyDiscreteContextToContinuousContext(discrete_context);
port.Calc(*(this->continuous_context_), out);
});
}
}

// Declare parameters.
for (int i = 0; i < continuous_context_->num_numeric_parameter_groups();
++i) {
this->DeclareNumericParameter(
continuous_context_->get_numeric_parameter(i));
}
for (int i = 0; i < continuous_context_->num_abstract_parameters(); ++i) {
this->DeclareAbstractParameter(
continuous_context_->get_abstract_parameter(i));
}

// Declare state.
this->DeclareDiscreteState(continuous_system_->num_continuous_states());
this->DeclarePeriodicDiscreteUpdateEvent(
time_period_, time_offset_, &DiscreteTimeSystem<T>::DiscreteUpdate);
}

void DiscreteUpdate(const Context<T>& discrete_context,
DiscreteValues<T>* out) const {
this->CopyDiscreteContextToContinuousContext(discrete_context);
this->integrator_->IntegrateWithMultipleStepsToTime(
continuous_context_->get_time() + time_period_);
out->get_mutable_vector().SetFromVector(
continuous_context_->get_continuous_state_vector().CopyToVector());
}

void CopyDiscreteContextToContinuousContext(
const Context<T>& discrete_context) const {
// Copy time.
continuous_context_->SetTime(discrete_context.get_time());
// Copy state.
continuous_context_->SetContinuousState(
discrete_context.get_discrete_state_vector().value());
// Copy parameters.
continuous_context_->get_mutable_parameters().SetFrom(
discrete_context.get_parameters());
// Copy accuracy.
continuous_context_->SetAccuracy(discrete_context.get_accuracy());
// Copy fixed input port values.
for (int i = 0; i < continuous_system_->num_input_ports(); ++i) {
continuous_context_->FixInputPort(
i, *(this->EvalAbstractInput(discrete_context, i)));
}
}

const std::unique_ptr<const System<T>> continuous_system_;
const std::unique_ptr<Context<T>> continuous_context_;
const double time_period_;
const double time_offset_;
const SimulatorConfig integrator_config_;
Simulator<T> simulator__;
IntegratorBase<T>* integrator_;
};

} // namespace

namespace scalar_conversion {
template <>
struct Traits<DiscreteTimeSystem> : public NonSymbolicTraits {};
} // namespace scalar_conversion

template <typename T>
std::unique_ptr<System<T>> DiscreteTimeApproximation(
const System<T>& system, double time_period, double time_offset,
const SimulatorConfig& integrator_config) {
// Check that the original system is continuous.
DRAKE_THROW_UNLESS(system.IsDifferentialEquationSystem());
// Check that the discrete time_period is greater than zero.
DRAKE_THROW_UNLESS(time_period > 0);

return std::make_unique<DiscreteTimeSystem<T>>(
system.Clone(), time_period, time_offset, integrator_config);
}

DRAKE_DEFINE_FUNCTION_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS(
(&DiscreteTimeApproximation<T>));

} // namespace systems
} // namespace drake
36 changes: 36 additions & 0 deletions systems/primitives/discrete_time_approximation.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,36 @@
#pragma once

#include <memory>

#include "drake/systems/analysis/simulator_config.h"
#include "drake/systems/framework/system.h"

namespace drake {
namespace systems {

/// Converts a general continuous-time system @f$ \dot{x} = f(t,x(t),u(t)) @f$
/// to a discrete-time system with zero-order hold on the input. The approximate
/// discrete-time dynamics is given by @f$ x[n+1] = f_d(n,x[n],u[n]) = x[n] +
/// \int_{t[n]}^{t[n+1]} f(t,x(t),u[n]) \, dt @f$, where the integration is
/// performed using numerical Integrators.
///
/// @param system The continuous-time System.
/// @param time_period The discrete time period.
/// @param time_period The discrete time offset.
/// @param integrator_config Use this parameter to choose a non-default
/// integrator or integrator parameters.
///
/// @returns A discrete-time System.
/// @pre @p system must be a continuous-time system. @p time_period must be
/// greater than zero.
///
/// @tparam_nonsymbolic_scalar
/// @ingroup primitive_systems
/// @pydrake_mkdoc_identifier{system}
template <typename T>
std::unique_ptr<System<T>> DiscreteTimeApproximation(
const System<T>& system, double time_period, double time_offset = 0.0,
const SimulatorConfig& integrator_config = SimulatorConfig());

} // namespace systems
} // namespace drake
Loading

0 comments on commit 6d44217

Please sign in to comment.