diff --git a/bindings/pybind11_ext/numpy_dtypes_user.h b/bindings/pybind11_ext/numpy_dtypes_user.h
index 983ce9335e8e..a12d60a3cd01 100644
--- a/bindings/pybind11_ext/numpy_dtypes_user.h
+++ b/bindings/pybind11_ext/numpy_dtypes_user.h
@@ -368,8 +368,8 @@ class dtype_user : public object {
   using Class = Class_;
   using DTypePyObject = detail::dtype_user_instance<Class>;
 
-  dtype_user(handle scope, const char* name) : cls_(none()) {
-    register_type(name);
+  dtype_user(handle scope, const char* name, const char* doc = "") : cls_(none()) {
+    register_type(name, doc);
     scope.attr(name) = self();
     auto& entry = detail::dtype_info::get_mutable_entry<Class>(true);
     entry.cls = self();
@@ -411,14 +411,14 @@ class dtype_user : public object {
   }
 
   template <typename ... Args>
-  dtype_user& def(detail::initimpl::constructor<Args...>&&) {
+  dtype_user& def(detail::initimpl::constructor<Args...>&&, const char* doc = "") {
     // See notes in `add_init`.
     // N.B. Do NOT use `Class*` as the argument, since that may incur recursion.
     add_init([](object py_self, Args... args) {
       // Old-style. No factories for now.
         Class* self = DTypePyObject::load_raw(py_self.ptr());
       new (self) Class(std::forward<Args>(args)...);
-    });
+    }, doc);
     return *this;
   }
 
@@ -547,7 +547,7 @@ class dtype_user : public object {
   }
 
   template <typename Func>
-  void add_init(Func&& f) {
+  void add_init(Func&& f, const char* doc) {
     // Do not construct this with the name `__init__` as `cpp_function` will
     // try to have this register the instance, and most likely segfault.
     this->def("_dtype_init", std::forward<Func>(f));
@@ -559,7 +559,7 @@ class dtype_user : public object {
           [init](handle self, args args, kwargs kwargs) {
             // Dispatch.
             init(self, *args, **kwargs);
-          }, is_method(self()));
+          }, is_method(self()), doc);
       self().attr("__init__") = func;
     }
   }
@@ -606,7 +606,7 @@ class dtype_user : public object {
     return 1;
   }
 
-  void register_type(const char* name) {
+  void register_type(const char* name, const char* doc) {
     // Ensure we initialize NumPy before accessing `PyGenericArrType_Type`.
     auto& api = detail::npy_api::get();
     // Loosely uses https://stackoverflow.com/a/12505371/7829525 as well.
@@ -624,6 +624,7 @@ class dtype_user : public object {
     ClassObject_Type.tp_basicsize = sizeof(DTypePyObject);
     ClassObject_Type.tp_getset = 0;
     ClassObject_Type.tp_flags = Py_TPFLAGS_DEFAULT | Py_TPFLAGS_HEAPTYPE;
+    ClassObject_Type.tp_doc = doc;
     if (PyType_Ready(&ClassObject_Type) != 0)
         pybind11_fail("dtype_user: Unable to initialize class");
     // TODO(eric.cousineau): Figure out how to catch recursions with
diff --git a/bindings/pydrake/_symbolic_extra.py b/bindings/pydrake/_symbolic_extra.py
index f4954aabba97..6a798e3d933e 100644
--- a/bindings/pydrake/_symbolic_extra.py
+++ b/bindings/pydrake/_symbolic_extra.py
@@ -2,6 +2,30 @@
 # rationale.
 import functools
 
+# Import aliases.
+# TODO(eric.cousineau): Deprecate, then remove these in lieu of `np.{func}`
+from pydrake.math import (
+    log,
+    abs,
+    exp,
+    pow,
+    sqrt,
+    sin,
+    cos,
+    tan,
+    asin,
+    acos,
+    atan,
+    atan2,
+    sinh,
+    cosh,
+    tanh,
+    min,
+    max,
+    ceil,
+    floor
+)
+
 
 def logical_and(*formulas):
     assert len(formulas) >= 1, "Must supply at least one operand"
diff --git a/bindings/pydrake/autodiffutils_py.cc b/bindings/pydrake/autodiffutils_py.cc
index 7004bcf0cecb..2e1b2b96f94c 100644
--- a/bindings/pydrake/autodiffutils_py.cc
+++ b/bindings/pydrake/autodiffutils_py.cc
@@ -29,103 +29,96 @@ PYBIND11_MODULE(autodiffutils, m) {
 
   py::dtype_user<AutoDiffXd> autodiff(m, "AutoDiffXd");
   DefImplicitConversionsFromNumericTypes(&autodiff);
-  autodiff
-    .def(py::init<double>())
-    .def(py::init<const double&, const Eigen::VectorXd&>())
-    // Downcasting must be explicit, to prevent inadvertent information loss.
-    .def_loop(py::dtype_method::explicit_conversion(
-        [](const AutoDiffXd& self) -> double { return self.value(); }))
-    // General methods.
-    .def("value", [](const AutoDiffXd& self) {
-      return self.value();
-    })
-    .def("derivatives", [](const AutoDiffXd& self) {
-      return self.derivatives();
-    })
-    .def("__str__", [](const AutoDiffXd& self) {
-      return py::str("AD{{{}, nderiv={}}}").format(
-          self.value(), self.derivatives().size());
-    })
-    .def("__repr__", [](const AutoDiffXd& self) {
-      return py::str("<AutoDiffXd {} nderiv={}>").format(
-          self.value(), self.derivatives().size());
-    })
-    // Arithmetic
-    .def_loop(-py::self)
-    .def_loop(py::self + py::self)
-    .def_loop(py::self + double())
-    .def_loop(double() + py::self)
-    .def_loop(py::self - py::self)
-    .def_loop(py::self - double())
-    .def_loop(double() - py::self)
-    .def_loop(py::self * py::self)
-    .def_loop(py::self * double())
-    .def_loop(double() * py::self)
-    .def_loop(py::self / py::self)
-    .def_loop(py::self / double())
-    .def_loop(double() / py::self)
-    // Logical comparison
-    .def_loop(py::self == py::self)
-    .def_loop(py::self == double())
-    .def_loop(py::self != py::self)
-    .def_loop(py::self != double())
-    // .def_loop(double() != py::self)
-    .def_loop(py::self < py::self)
-    .def_loop(py::self < double())
-    .def_loop(double() < py::self)
-    .def_loop(py::self <= py::self)
-    .def_loop(py::self <= double())
-    .def_loop(double() <= py::self)
-    .def_loop(py::self > py::self)
-    .def_loop(py::self > double())
-    .def_loop(double() > py::self)
-    .def_loop(py::self >= py::self)
-    .def_loop(py::self >= double())
-    .def_loop(double() >= py::self)
-    // Dot-product
-    .def_loop(py::dtype_method::dot())
-    // N.B. Certain versions of NumPy require `square` be specifically
-    // defined.
-    .def_loop("square", [](const AutoDiffXd& self) {
-      return self * self;
-    })
-    .def_loop("isfinite", [](const AutoDiffXd& self) {
-      return isfinite(self.value());
-    });
+  autodiff  // BR
+      .def(py::init<double>())
+      .def(py::init<const double&, const Eigen::VectorXd&>())
+      // Downcasting must be explicit, to prevent inadvertent information loss.
+      .def_loop(py::dtype_method::explicit_conversion(
+          [](const AutoDiffXd& self) -> double { return self.value(); }))
+      // General methods.
+      .def("value", [](const AutoDiffXd& self) { return self.value(); })
+      .def("derivatives",
+           [](const AutoDiffXd& self) { return self.derivatives(); })
+      .def("__str__",
+           [](const AutoDiffXd& self) {
+             return py::str("AD{{{}, nderiv={}}}")
+                 .format(self.value(), self.derivatives().size());
+           })
+      .def("__repr__",
+           [](const AutoDiffXd& self) {
+             return py::str("<AutoDiffXd {} nderiv={}>")
+                 .format(self.value(), self.derivatives().size());
+           })
+      // Arithmetic
+      .def_loop(-py::self)
+      .def_loop(py::self + py::self)
+      .def_loop(py::self + double())
+      .def_loop(double() + py::self)
+      .def_loop(py::self - py::self)
+      .def_loop(py::self - double())
+      .def_loop(double() - py::self)
+      .def_loop(py::self * py::self)
+      .def_loop(py::self * double())
+      .def_loop(double() * py::self)
+      .def_loop(py::self / py::self)
+      .def_loop(py::self / double())
+      .def_loop(double() / py::self)
+      // Logical comparison
+      .def_loop(py::self == py::self)
+      .def_loop(py::self == double())
+      .def_loop(py::self != py::self)
+      .def_loop(py::self != double())
+      // .def_loop(double() != py::self)
+      .def_loop(py::self < py::self)
+      .def_loop(py::self < double())
+      .def_loop(double() < py::self)
+      .def_loop(py::self <= py::self)
+      .def_loop(py::self <= double())
+      .def_loop(double() <= py::self)
+      .def_loop(py::self > py::self)
+      .def_loop(py::self > double())
+      .def_loop(double() > py::self)
+      .def_loop(py::self >= py::self)
+      .def_loop(py::self >= double())
+      .def_loop(double() >= py::self)
+      // Dot-product
+      .def_loop(py::dtype_method::dot())
+      // N.B. Certain versions of NumPy require `square` be specifically
+      // defined.
+      .def_loop("square", [](const AutoDiffXd& self) { return self * self; })
+      .def_loop("isfinite",
+                [](const AutoDiffXd& self) { return isfinite(self.value()); });
 
   // Add overloads for `math` functions.
   auto math = py::module::import("pydrake.math");
   UfuncMirrorDef<decltype(autodiff)>(&autodiff, math)
-    .def_loop("__pow__", "pow",
-         [](const AutoDiffXd& base, int exponent) {
-           return pow(base, exponent);
-         })
-    .def_loop("__abs__", "abs", [](const AutoDiffXd& x) { return abs(x); })
-    .def_loop("log", [](const AutoDiffXd& x) { return log(x); })
-    .def_loop("exp", [](const AutoDiffXd& x) { return exp(x); })
-    .def_loop("sqrt", [](const AutoDiffXd& x) { return sqrt(x); })
-    .def_loop("sin", [](const AutoDiffXd& x) { return sin(x); })
-    .def_loop("cos", [](const AutoDiffXd& x) { return cos(x); })
-    .def_loop("tan", [](const AutoDiffXd& x) { return tan(x); })
-    .def_loop("arcsin", "asin", [](const AutoDiffXd& x) { return asin(x); })
-    .def_loop("arccos", "acos", [](const AutoDiffXd& x) { return acos(x); })
-    .def_loop("arctan2", "atan2",
-        [](const AutoDiffXd& y, const AutoDiffXd& x) {
-            return atan2(y, x);
-        })
-    .def_loop("sinh", [](const AutoDiffXd& x) { return sinh(x); })
-    .def_loop("cosh", [](const AutoDiffXd& x) { return cosh(x); })
-    .def_loop("tanh", [](const AutoDiffXd& x) { return tanh(x); })
-    .def_loop("fmin", "min",
-        [](const AutoDiffXd& x, const AutoDiffXd& y) {
-            return min(x, y);
-        })
-    .def_loop("fmax", "max",
-        [](const AutoDiffXd& x, const AutoDiffXd& y) {
-            return max(x, y);
-        })
-    .def_loop("ceil", [](const AutoDiffXd& x) { return ceil(x); })
-    .def_loop("floor", [](const AutoDiffXd& x) { return floor(x); });
+      .def_loop("__pow__", "pow",
+                [](const AutoDiffXd& base, int exponent) {
+                  return pow(base, exponent);
+                })
+      .def_loop("__abs__", "abs", [](const AutoDiffXd& x) { return abs(x); })
+      .def_loop("log", [](const AutoDiffXd& x) { return log(x); })
+      .def_loop("exp", [](const AutoDiffXd& x) { return exp(x); })
+      .def_loop("sqrt", [](const AutoDiffXd& x) { return sqrt(x); })
+      .def_loop("sin", [](const AutoDiffXd& x) { return sin(x); })
+      .def_loop("cos", [](const AutoDiffXd& x) { return cos(x); })
+      .def_loop("tan", [](const AutoDiffXd& x) { return tan(x); })
+      .def_loop("arcsin", "asin", [](const AutoDiffXd& x) { return asin(x); })
+      .def_loop("arccos", "acos", [](const AutoDiffXd& x) { return acos(x); })
+      .def_loop(
+          "arctan2", "atan2",
+          [](const AutoDiffXd& y, const AutoDiffXd& x) { return atan2(y, x); })
+      .def_loop("sinh", [](const AutoDiffXd& x) { return sinh(x); })
+      .def_loop("cosh", [](const AutoDiffXd& x) { return cosh(x); })
+      .def_loop("tanh", [](const AutoDiffXd& x) { return tanh(x); })
+      .def_loop(
+          "fmin", "min",
+          [](const AutoDiffXd& x, const AutoDiffXd& y) { return min(x, y); })
+      .def_loop(
+          "fmax", "max",
+          [](const AutoDiffXd& x, const AutoDiffXd& y) { return max(x, y); })
+      .def_loop("ceil", [](const AutoDiffXd& x) { return ceil(x); })
+      .def_loop("floor", [](const AutoDiffXd& x) { return floor(x); });
 }
 
 }  // namespace pydrake
diff --git a/bindings/pydrake/pydrake_pybind.h b/bindings/pydrake/pydrake_pybind.h
index d2c21146597c..e0cb7287bdbb 100644
--- a/bindings/pydrake/pydrake_pybind.h
+++ b/bindings/pydrake/pydrake_pybind.h
@@ -276,7 +276,8 @@ This is a brief recipe for debugging with GDB
         source_dir=${workspace}/bazel-${name}
         source /tmp/env.sh
         cd ${target_bin_path}.runfiles/${name}
-        gdb --directory ${source_dir} --args python ${target_bin_path} --trace=user -f
+        gdb --directory ${source_dir} --args python ${target_bin_path}
+--trace=user -f
     )
 
 This allows you to use GDB from the terminal, while being able to inspect the
@@ -333,7 +334,8 @@ template <typename PyClass>
 void DefCopyAndDeepCopy(PyClass* ppy_class) {
   using Class = typename PyClass::type;
   PyClass& py_class = *ppy_class;
-  py_class.def("__copy__", [](const Class* self) { return Class{*self}; })
+  py_class  // BR
+      .def("__copy__", [](const Class* self) { return Class{*self}; })
       .def("__deepcopy__",
            [](const Class* self, py::dict /* memo */) { return Class{*self}; });
 }
diff --git a/bindings/pydrake/symbolic_py.cc b/bindings/pydrake/symbolic_py.cc
index 06ceb2a28944..b67a60467815 100644
--- a/bindings/pydrake/symbolic_py.cc
+++ b/bindings/pydrake/symbolic_py.cc
@@ -46,7 +46,7 @@ PYBIND11_MODULE(symbolic, m) {
   py::dtype_user<Expression> expr_cls(m, "Expression", doc.Expression.doc);
   py::dtype_user<Formula> formula_cls(m, "Formula", doc.Formula.doc);
 
-  var_cls
+  var_cls  // BR
       .def(py::init<const string&>(), doc.Variable.ctor.doc)
       .def("get_id", &Variable::get_id, doc.Variable.get_id.doc)
       .def("__str__", &Variable::to_string, doc.Variable.to_string.doc)
@@ -76,20 +76,18 @@ PYBIND11_MODULE(symbolic, m) {
       // `dot` for NumPy since it expects a closed operation.
       // TODO(eric.cousineau): See if `dot` can be defined at some point.
       // Pow.
+      .def_loop("__pow__", [](const Variable& self,
+                              double other) { return pow(self, other); })
       .def_loop("__pow__",
-           [](const Variable& self, double other) { return pow(self, other); })
+                [](const Variable& self, const Variable& other) {
+                  return pow(self, other);
+                })
       .def_loop("__pow__",
-           [](const Variable& self, const Variable& other) {
-             return pow(self, other);
-           })
-      .def_loop("__pow__",
-           [](const Variable& self, const Expression& other) {
-             return pow(self, other);
-           })
+                [](const Variable& self, const Expression& other) {
+                  return pow(self, other);
+                })
       // See comment about `np.square` in AutoDiff<> bindings.
-      .def_loop("square", [](const Variable& self) {
-        return self * self;
-      })
+      .def_loop("square", [](const Variable& self) { return self * self; })
       // We add `EqualTo` instead of `equal_to` to maintain consistency among
       // symbolic classes (Variable, Expression, Formula, Polynomial) on Python
       // side. This enables us to achieve polymorphism via ducktyping in Python.
@@ -128,7 +126,7 @@ PYBIND11_MODULE(symbolic, m) {
       .def_loop(py::self != Expression())
       .def_loop(py::self != py::self)
       .def_loop(py::self != double());
-  DefCopyAndDeepCopy(&var_cls);
+  DefCopyAndDeepCopy(&var_cls.cls());
 
   // TODO(m-chaturvedi) Add Pybind11 documentation for operator overloads, etc.
   py::class_<Variables>(m, "Variables", doc.Variables.doc)
@@ -191,7 +189,7 @@ PYBIND11_MODULE(symbolic, m) {
   });
 
   DefImplicitConversionsFromNumericTypes(&expr_cls);
-  expr_cls
+  expr_cls  // BR
       .def(py::init<>(), doc.Expression.ctor.doc_3)
       .def(py::init<double>(), doc.Expression.ctor.doc_4)
       .def(py::init<const Variable&>(), doc.Expression.ctor.doc_5)
@@ -304,9 +302,7 @@ PYBIND11_MODULE(symbolic, m) {
       .def_loop(py::self != Variable())
       .def_loop(py::self != double())
       // See comment about `np.square` in AutoDiff<> bindings.
-      .def_loop("square", [](const Expression& self) {
-        return self * self;
-      })
+      .def_loop("square", [](const Expression& self) { return self * self; })
       .def("Differentiate", &Expression::Differentiate,
            doc.Expression.Differentiate.doc)
       .def("Jacobian", &Expression::Jacobian);
@@ -325,17 +321,22 @@ PYBIND11_MODULE(symbolic, m) {
       // TODO(eric.cousineau): Figure out how to consolidate with the below
       // methods.
       // Pow.
-      .def_loop("__pow__", "pow", [](const Expression& self,
-                         const double other) { return pow(self, other); })
-      .def_loop("__pow__", "pow", [](const Expression& self,
-                         const Variable& other) { return pow(self, other); })
-      .def_loop("__pow__", "pow", [](const Expression& self,
-                         const Expression& other) { return pow(self, other); })
+      .def_loop("__pow__", "pow",
+                [](const Expression& self, const double other) {
+                  return pow(self, other);
+                })
+      .def_loop("__pow__", "pow",
+                [](const Expression& self, const Variable& other) {
+                  return pow(self, other);
+                })
+      .def_loop("__pow__", "pow",
+                [](const Expression& self, const Expression& other) {
+                  return pow(self, other);
+                })
       .def_loop("log", &symbolic::log)
       .def_loop("__abs__", "abs", &symbolic::abs)
       .def_loop("exp", &symbolic::exp)
       .def_loop("sqrt", &symbolic::sqrt)
-      // TODO(eric.cousineau): Move `__pow__` here.
       .def_loop("sin", &symbolic::sin)
       .def_loop("cos", &symbolic::cos)
       .def_loop("tan", &symbolic::tan)
@@ -353,32 +354,6 @@ PYBIND11_MODULE(symbolic, m) {
   MirrorDef<decltype(expr_cls), py::module>(&expr_cls, &math)
       .def("atan", &symbolic::atan);
 
-  // Import aliases.
-  // TODO(eric.cousineau): Deprecate, then remove these in lieu of `np.{func}`
-  py::exec(R"""(
-from pydrake.math import (
-    log,
-    abs,
-    exp,
-    pow,
-    sqrt,
-    sin,
-    cos,
-    tan,
-    asin,
-    acos,
-    atan,
-    atan2,
-    sinh,
-    cosh,
-    tanh,
-    min,
-    max,
-    ceil,
-    floor
-)
-)""");
-
   m.def("if_then_else", [](bool cond, double true_value, double false_value) {
     return cond ? true_value : false_value;
   });
@@ -426,8 +401,10 @@ from pydrake.math import (
            [](const Formula& self) {
              return fmt::format("<Formula \"{}\">", self.to_string());
            })
-      .def_loop("__eq__", [](const Formula& self,
-                        const Formula& other) { return self.EqualTo(other); })
+      .def_loop("__eq__",
+                [](const Formula& self, const Formula& other) {
+                  return self.EqualTo(other);
+                })
       .def("__ne__", [](const Formula& self,
                         const Formula& other) { return !self.EqualTo(other); })
       .def("__hash__",
diff --git a/bindings/pydrake/systems/framework_py_semantics.cc b/bindings/pydrake/systems/framework_py_semantics.cc
index 4976d77aadef..7b88f7ef207e 100644
--- a/bindings/pydrake/systems/framework_py_semantics.cc
+++ b/bindings/pydrake/systems/framework_py_semantics.cc
@@ -347,6 +347,7 @@ void DefineFrameworkPySemantics(py::module m) {
              // Keep alive, ownership: `value` keeps `self` alive.
              py::keep_alive<2, 1>(), py::arg("abstract_params"),
              doc.Parameters.set_abstract_parameters.doc)
+        .def("CopyFrom", &Parameters<T>::CopyFrom, doc.Parameters.CopyFrom.doc)
         .def("SetFrom", &Parameters<T>::SetFrom, doc.Parameters.SetFrom.doc);
 
     // State.
diff --git a/bindings/pydrake/systems/framework_py_systems.cc b/bindings/pydrake/systems/framework_py_systems.cc
index b29b533acdad..2ef58fdf9b26 100644
--- a/bindings/pydrake/systems/framework_py_systems.cc
+++ b/bindings/pydrake/systems/framework_py_systems.cc
@@ -497,10 +497,11 @@ struct Impl {
                return new PyVectorSystem(inputs, outputs);
              }),
              doc.VectorSystem.ctor.doc)
-        .def(py::init(
-            [](int inputs, int outputs, SystemScalarConverter converter) {
-              return new PyVectorSystem(std::move(converter), inputs, outputs);
-            }), py::arg("inputs"), py::arg("outputs"), py::arg("converter"));
+        .def(py::init([](int inputs, int outputs,
+                         SystemScalarConverter converter) {
+               return new PyVectorSystem(std::move(converter), inputs, outputs);
+             }),
+             py::arg("inputs"), py::arg("outputs"), py::arg("converter"));
     // TODO(eric.cousineau): Bind virtual methods once we provide a function
     // wrapper to convert `Map<Derived>*` arguments.
     // N.B. This could be mitigated by using `EigenPtr` in public interfaces in
diff --git a/bindings/pydrake/systems/framework_py_values.cc b/bindings/pydrake/systems/framework_py_values.cc
index a88d0cf85f0c..df90243c601a 100644
--- a/bindings/pydrake/systems/framework_py_values.cc
+++ b/bindings/pydrake/systems/framework_py_values.cc
@@ -48,15 +48,15 @@ void DefineFrameworkPyValues(py::module m) {
             m, "BasicVector", GetPyParam<T>(), doc.BasicVector.doc);
     DefClone(&basic_vector);
     basic_vector
-      // N.B. Place `init<VectorX<T>>` `init<int>` so that we do not implicitly
-      // convert scalar-size `np.array` objects to `int` (since this is normally
-      // permitted).
-      // N.B. Also ensure that we use `greedy_arg` to prevent ambiguous
-      // overloads when using scalars vs. lists vs. numpy arrays. See
-      // `greedy_arg` for more information.
-      .def(py::init([](greedy_arg<VectorX<T>> in) {
-        return new BasicVector<T>(*in);
-      }), doc.BasicVector.ctor.doc_5)
+        // N.B. Place `init<VectorX<T>>` `init<int>` so that we do not
+        // implicitly convert scalar-size `np.array` objects to `int` (since
+        // this is normally permitted). N.B. Also ensure that we use
+        // `greedy_arg` to prevent ambiguous overloads when using scalars vs.
+        // lists vs. numpy arrays. See `greedy_arg` for more information.
+        .def(py::init([](greedy_arg<VectorX<T>> in) {
+               return new BasicVector<T>(*in);
+             }),
+             doc.BasicVector.ctor.doc_5)
         .def(py::init<int>(), py::arg("size"), doc.BasicVector.ctor.doc_4)
         .def("get_value",
              [](const BasicVector<T>* self) -> Eigen::Ref<const VectorX<T>> {
diff --git a/bindings/pydrake/systems/test/custom_test.py b/bindings/pydrake/systems/test/custom_test.py
index a2f10c93299c..d47c0a567adb 100644
--- a/bindings/pydrake/systems/test/custom_test.py
+++ b/bindings/pydrake/systems/test/custom_test.py
@@ -14,12 +14,12 @@
 from pydrake.systems.framework import (
     AbstractValue,
     BasicVector, BasicVector_,
-    Context,
+    Context, Context_,
     DiagramBuilder, DiagramBuilder_,
     kUseDefaultName,
     LeafSystem, LeafSystem_,
     PortDataType,
-    PublishEvent,
+    PublishEvent, PublishEvent_,
     TriggerType,
     VectorSystem, VectorSystem_,
     )
@@ -187,7 +187,7 @@ def __init__(self):
                 self._DeclarePeriodicDiscreteUpdate(
                     period_sec=0.1, offset_sec=0.)
                 self._DeclareInitializationEvent(
-                    event=PublishEvent(
+                    event=PublishEvent_[T](
                         trigger_type=TriggerType.kInitialization,
                         callback=self._on_initialize))
                 self._DeclareContinuousState(2)
@@ -234,8 +234,8 @@ def _DoCalcDiscreteVariableUpdates(
                 self.called_discrete = True
 
             def _on_initialize(self, context, event):
-                test.assertIsInstance(context, Context)
-                test.assertIsInstance(event, PublishEvent)
+                test.assertIsInstance(context, Context_[T])
+                test.assertIsInstance(event, PublishEvent_[T])
                 test.assertFalse(self.called_initialize)
                 self.called_initialize = True
 
diff --git a/bindings/pydrake/systems/test/test_util_py.cc b/bindings/pydrake/systems/test/test_util_py.cc
index b3d25e425eb5..74b3dfbc9afe 100644
--- a/bindings/pydrake/systems/test/test_util_py.cc
+++ b/bindings/pydrake/systems/test/test_util_py.cc
@@ -14,8 +14,10 @@ using std::unique_ptr;
 namespace drake {
 
 using systems::BasicVector;
+using systems::Context;
 using systems::LeafSystem;
 using systems::Simulator;
+using systems::System;
 
 namespace pydrake {
 namespace {
@@ -79,6 +81,23 @@ class MyVector2 : public BasicVector<T> {
   }
 };
 
+// Stolen from `Simulator::Initialize`.
+// TODO(eric.cousineau): Somehow leverage `Simulator` bits, factoring that
+// functionality out.
+template <typename T>
+void DispatchInitializationEvents(const System<T>& system,
+                                  const Context<T>& context) {
+  // Process all the initialization events.
+  auto init_events = system.AllocateCompositeEventCollection();
+  system.GetInitializationEvents(context, init_events.get());
+  DRAKE_DEMAND(!init_events->get_unrestricted_update_events().HasEvents());
+  DRAKE_DEMAND(!init_events->get_discrete_update_events().HasEvents());
+  auto& pub_events = init_events->get_publish_events();
+  if (pub_events.HasEvents()) {
+    system.Publish(context, pub_events);
+  }
+}
+
 }  // namespace
 
 PYBIND11_MODULE(test_util, m) {
@@ -119,9 +138,12 @@ PYBIND11_MODULE(test_util, m) {
     };
 
     m.def("call_leaf_system_overrides", [clone_vector](
-        const LeafSystem<T>& system) {
+                                            const LeafSystem<T>& system) {
       py::dict results;
       auto context = system.AllocateContext();
+      // Check initialization events.
+      DispatchInitializationEvents(system, *context);
+
       {
         // Call `Publish` to test `DoPublish`.
         auto events =
@@ -149,113 +171,40 @@ PYBIND11_MODULE(test_util, m) {
         // If there is an abstract / unrestricted update, this assumes that
         // `dt_discrete < dt_abstract`.
         systems::LeafCompositeEventCollection<T> events;
-        results["discrete_next_t"] = system.CalcNextUpdateTime(
-            *context, &events);
+        results["discrete_next_t"] =
+            system.CalcNextUpdateTime(*context, &events);
       }
       return results;
     });
 
-    m.def("call_vector_system_overrides", [clone_vector](
-        const VectorSystem<T>& system, Context<T>* context,
-        bool is_discrete, double dt) {
-      // While this is not convention, update state first to ensure that our
-      // output incorporates it correctly, for testing purposes.
-      // TODO(eric.cousineau): Add (Continuous|Discrete)State::Clone().
-      if (is_discrete) {
-        auto& state = context->get_mutable_discrete_state();
-        DiscreteValues<T> state_copy(
-            clone_vector(state.get_vector()));
-        system.CalcDiscreteVariableUpdates(
-            *context, &state_copy);
-      state.CopyFrom(state_copy);
-      } else {
-        auto& state = context->get_mutable_continuous_state();
-        ContinuousState<T> state_dot(
-            clone_vector(state.get_vector()),
-            state.get_generalized_position().size(),
-            state.get_generalized_velocity().size(),
-            state.get_misc_continuous_state().size());
-        system.CalcTimeDerivatives(*context, &state_dot);
-        state.SetFromVector(
-            state.CopyToVector() + dt * state_dot.CopyToVector());
-      }
-      // Calculate output.
-      auto output = system.AllocateOutput();
-      system.CalcOutput(*context, output.get());
-      return output;
-    });
+    m.def("call_vector_system_overrides",
+          [clone_vector](const VectorSystem<T>& system, Context<T>* context,
+                         bool is_discrete, double dt) {
+            // While this is not convention, update state first to ensure that
+            // our output incorporates it correctly, for testing purposes.
+            // TODO(eric.cousineau): Add (Continuous|Discrete)State::Clone().
+            if (is_discrete) {
+              auto& state = context->get_mutable_discrete_state();
+              DiscreteValues<T> state_copy(clone_vector(state.get_vector()));
+              system.CalcDiscreteVariableUpdates(*context, &state_copy);
+              state.CopyFrom(state_copy);
+            } else {
+              auto& state = context->get_mutable_continuous_state();
+              ContinuousState<T> state_dot(
+                  clone_vector(state.get_vector()),
+                  state.get_generalized_position().size(),
+                  state.get_generalized_velocity().size(),
+                  state.get_misc_continuous_state().size());
+              system.CalcTimeDerivatives(*context, &state_dot);
+              state.SetFromVector(state.CopyToVector() +
+                                  dt * state_dot.CopyToVector());
+            }
+            // Calculate output.
+            auto output = system.AllocateOutput();
+            system.CalcOutput(*context, output.get());
+            return output;
+          });
   };
-
-  m.def("call_leaf_system_overrides", [clone_vector](
-                                          const LeafSystem<T>& system) {
-    py::dict results;
-    auto context = system.AllocateContext();
-    {
-      // Leverage simulator to call initialization events.
-      Simulator<T> simulator(system);
-      // Do not publish at initialization because we want to track publishes
-      // from only events of trigger type `kInitialization`.
-      simulator.set_publish_at_initialization(false);
-      simulator.Initialize();
-    }
-    {
-      // Call `Publish` to test `DoPublish`.
-      auto events =
-          LeafEventCollection<PublishEvent<T>>::MakeForcedEventCollection();
-      system.Publish(*context, *events);
-    }
-    {
-      // Call `HasDirectFeedthrough` to test `DoHasDirectFeedthrough`.
-      results["has_direct_feedthrough"] = system.HasDirectFeedthrough(0, 0);
-    }
-    {
-      // Call `CalcTimeDerivatives` to test `DoCalcTimeDerivatives`
-      auto& state = context->get_mutable_continuous_state();
-      ContinuousState<T> state_copy(clone_vector(state.get_vector()));
-      system.CalcTimeDerivatives(*context, &state_copy);
-    }
-    {
-      // Call `CalcDiscreteVariableUpdates` to test
-      // `DoCalcDiscreteVariableUpdates`.
-      auto& state = context->get_mutable_discrete_state();
-      DiscreteValues<T> state_copy(clone_vector(state.get_vector()));
-      system.CalcDiscreteVariableUpdates(*context, &state_copy);
-
-      // From t=0, return next update time for testing discrete time.
-      // If there is an abstract / unrestricted update, this assumes that
-      // `dt_discrete < dt_abstract`.
-      systems::LeafCompositeEventCollection<double> events;
-      results["discrete_next_t"] = system.CalcNextUpdateTime(*context, &events);
-    }
-    return results;
-  });
-
-  m.def("call_vector_system_overrides", [clone_vector](
-                                            const VectorSystem<T>& system,
-                                            Context<T>* context,
-                                            bool is_discrete, double dt) {
-    // While this is not convention, update state first to ensure that our
-    // output incorporates it correctly, for testing purposes.
-    // TODO(eric.cousineau): Add (Continuous|Discrete)State::Clone().
-    if (is_discrete) {
-      auto& state = context->get_mutable_discrete_state();
-      DiscreteValues<T> state_copy(clone_vector(state.get_vector()));
-      system.CalcDiscreteVariableUpdates(*context, &state_copy);
-      state.SetFrom(state_copy);
-    } else {
-      auto& state = context->get_mutable_continuous_state();
-      ContinuousState<T> state_dot(clone_vector(state.get_vector()),
-                                   state.get_generalized_position().size(),
-                                   state.get_generalized_velocity().size(),
-                                   state.get_misc_continuous_state().size());
-      system.CalcTimeDerivatives(*context, &state_dot);
-      state.SetFromVector(state.CopyToVector() + dt * state_dot.CopyToVector());
-    }
-    // Calculate output.
-    auto output = system.AllocateOutput();
-    system.CalcOutput(*context, output.get());
-    return output;
-  });
   type_visit(bind_common_scalar_types, pysystems::CommonScalarPack{});
 }
 
diff --git a/bindings/pydrake/test/autodiffutils_test_util_py.cc b/bindings/pydrake/test/autodiffutils_test_util_py.cc
index 89bfc654eed6..db311411298e 100644
--- a/bindings/pydrake/test/autodiffutils_test_util_py.cc
+++ b/bindings/pydrake/test/autodiffutils_test_util_py.cc
@@ -13,11 +13,11 @@ class AutoDiffContainer {
  public:
   AutoDiffContainer() {
     value_.resize(1, 2);
-    value_ <<
-        AutoDiffXd(10, (VectorXd(2) << 1, 0).finished()),
+    value_ << AutoDiffXd(10, (VectorXd(2) << 1, 0).finished()),
         AutoDiffXd(100, (VectorXd(2) << 0, 1).finished());
   }
   MatrixX<AutoDiffXd>& value() { return value_; }
+
  private:
   MatrixX<AutoDiffXd> value_;
 };
@@ -28,22 +28,18 @@ PYBIND11_MODULE(autodiffutils_test_util, m) {
   py::module::import("pydrake.autodiffutils");
 
   // Implicit argument conversion.
-  m.def("autodiff_scalar_pass_through", [](const AutoDiffXd& value) {
-    return value;
-  });
-  m.def("autodiff_vector_pass_through", [](const VectorX<AutoDiffXd>& value) {
-    return value;
-  });
+  m.def("autodiff_scalar_pass_through",
+        [](const AutoDiffXd& value) { return value; });
+  m.def("autodiff_vector_pass_through",
+        [](const VectorX<AutoDiffXd>& value) { return value; });
 
   // Reference semantics for AutoDiff.
   py::class_<AutoDiffContainer>(m, "AutoDiffContainer")
       .def(py::init())
-      .def("value", &AutoDiffContainer::value,
-           py_reference_internal);
+      .def("value", &AutoDiffContainer::value, py_reference_internal);
 
-  m.def("autodiff_increment", [](Eigen::Ref<MatrixX<AutoDiffXd>> value) {
-      value.array() += 1;
-  });
+  m.def("autodiff_increment",
+        [](Eigen::Ref<MatrixX<AutoDiffXd>> value) { value.array() += 1; });
 }
 
 }  // namespace pydrake
diff --git a/bindings/pydrake/test/symbolic_test.py b/bindings/pydrake/test/symbolic_test.py
index 27ed4681f34b..43ea28689ecb 100644
--- a/bindings/pydrake/test/symbolic_test.py
+++ b/bindings/pydrake/test/symbolic_test.py
@@ -366,7 +366,7 @@ def _check_algebra(self, algebra):
         algebra.check_value((1 - e_xv), "(1 - x)")
 
         # - In place.
-        e = copy(e_xv)  # N.B. No array construnction needed
+        e = copy.copy(e_xv)  # N.B. No array construnction needed
         e -= e_yv
         algebra.check_value(e, (x - y))
         e -= zv
@@ -676,8 +676,8 @@ def test_substitute_with_dict(self):
         self.assertEqualStructure(e.Substitute(env), x + y + 5)
 
     def test_copy(self):
-        self._check_scalar(copy.copy(e_x), e_x)
-        self._check_scalar(copy.deepcopy(e_x), e_x)
+        self.assertEqualStructure(copy.copy(e_x), e_x)
+        self.assertEqualStructure(copy.deepcopy(e_x), e_x)
 
     # See `math_overloads_test` for more comprehensive checks on math
     # functions.
diff --git a/tools/lint/clang_format_lint.py b/tools/lint/clang_format_lint.py
index f30345db56fa..90c42a3124e2 100644
--- a/tools/lint/clang_format_lint.py
+++ b/tools/lint/clang_format_lint.py
@@ -33,7 +33,8 @@ def _check_clang_format_idempotence(filename):
     if not changes:
         return 0
     print("ERROR: {} needs clang-format".format(filename))
-    print("note: fix via {} -style=file -i {}".format(clang_format, filename))
+    print("note: fix via\n  {} -style=file -i {}".format(
+        clang_format, filename))
     return 1
 
 
diff --git a/tools/skylark/pybind.bzl b/tools/skylark/pybind.bzl
index e0088aa646d9..8edb159bc2dd 100644
--- a/tools/skylark/pybind.bzl
+++ b/tools/skylark/pybind.bzl
@@ -73,6 +73,10 @@ def pybind_py_library(
         visibility = visibility,
     )
 
+    # Ensure that we're always using our custom NumPy version.
+    if "@numpy_py" not in py_deps:
+        py_deps = py_deps + ["@numpy_py"]
+
     # Add Python library.
     py_library_rule(
         name = py_name,