Hi everyone,
I am trying to solve an optimal control problem with an ode as a state equation. when simply solving the ode, I get the solution over all the vertices of the domain; but when I try to solve the optimal control problem, I am getting an error. Below is the code
import numpy as np
from dolfin import*
from dolfin_adjoint import*
import matplotlib.pyplot as plt
from collections import OrderedDict
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1) # Linear function space
time = Constant(0.0)
time = float(time)
final_time = 2.0
dt = Constant(0.1)
ctrls = OrderedDict()
tt = []
while time <= final_time:
ctrls[time] = Function(V)
tt.append(time)
time = round( (time + float(dt)), 1 )
y_desired = Function(V, name="desired state")
y = TrialFunction(V)
v = TestFunction(V)
t = Constant(0.0)
form = 2*t*ctrls[float(t)]*v*dP # variational form
y = project( Constant(2.0), V )
scheme = RK4(form, y, t)
ode_solver = PointIntegralSolver(scheme)
J = 0
num_step = int(final_time/float(dt))
for i in range(num_step):
t1 = round( (float(t) + float(dt)), 1 )
ode_solver.step(dt)
J += 0.5*float(dt)*assemble( (y-y_desired)**2*dP ) # continuous time integral discretized an we obtain a discrete sum
t.assign(t1)
m = [Control(ctrls[i]) for i in tt]
reduced_functional = ReducedFunctional(J, m)
m_opt = minimize(reduced_functional, method='L-BFGS-B')
Here is the error I am getting:
Not an UFL type: <class 'ufl.form.Form'>
---------------------------------------------------------------------------
UFLException Traceback (most recent call last)
<ipython-input-15-c180f0adc09d> in <module>
----> 1 m_opt = minimize(reduced_functional, method='L-BFGS-B')
~/.local/lib/python3.9/site-packages/pyadjoint/tape.py in wrapper(*args, **kwargs)
107 def wrapper(*args, **kwargs):
108 with stop_annotating():
--> 109 return function(*args, **kwargs)
110
111 return wrapper
~/.local/lib/python3.9/site-packages/pyadjoint/optimization/optimization.py in minimize(rf, method, scale, **kwargs)
249 kwargs["method"] = method
250
--> 251 opt = algorithm(rf_np, **kwargs)
252
253 if len(opt) == 1:
~/.local/lib/python3.9/site-packages/pyadjoint/optimization/optimization.py in minimize_scipy_generic(rf_np, method, bounds, **kwargs)
132 res = scipy_minimize(J, m_global, method=method, bounds=bounds, **kwargs)
133 else:
--> 134 res = scipy_minimize(J, m_global, method=method, **kwargs)
135
136 m = rf_np.set_controls(np.array(res["x"]))
/usr/lib/python3/dist-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
617 **options)
618 elif meth == 'l-bfgs-b':
--> 619 return _minimize_lbfgsb(fun, x0, args, jac, bounds,
620 callback=callback, **options)
621 elif meth == 'tnc':
/usr/lib/python3/dist-packages/scipy/optimize/lbfgsb.py in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, **unknown_options)
304 iprint = disp
305
--> 306 sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,
307 bounds=new_bounds,
308 finite_diff_rel_step=finite_diff_rel_step)
/usr/lib/python3/dist-packages/scipy/optimize/optimize.py in _prepare_scalar_function(fun, x0, jac, args, bounds, epsilon, finite_diff_rel_step, hess)
259 # ScalarFunction caches. Reuse of fun(x) during grad
260 # calculation reduces overall function evaluations.
--> 261 sf = ScalarFunction(fun, x0, args, grad, hess,
262 finite_diff_rel_step, bounds, epsilon=epsilon)
263
/usr/lib/python3/dist-packages/scipy/optimize/_differentiable_functions.py in __init__(self, fun, x0, args, grad, hess, finite_diff_rel_step, finite_diff_bounds, epsilon)
134
135 self._update_fun_impl = update_fun
--> 136 self._update_fun()
137
138 # Gradient evaluation
/usr/lib/python3/dist-packages/scipy/optimize/_differentiable_functions.py in _update_fun(self)
224 def _update_fun(self):
225 if not self.f_updated:
--> 226 self._update_fun_impl()
227 self.f_updated = True
228
/usr/lib/python3/dist-packages/scipy/optimize/_differentiable_functions.py in update_fun()
131
132 def update_fun():
--> 133 self.f = fun_wrapped(self.x)
134
135 self._update_fun_impl = update_fun
/usr/lib/python3/dist-packages/scipy/optimize/_differentiable_functions.py in fun_wrapped(x)
128 def fun_wrapped(x):
129 self.nfev += 1
--> 130 return fun(x, *args)
131
132 def update_fun():
~/.local/lib/python3.9/site-packages/pyadjoint/reduced_functional_numpy.py in __call__(self, m_array)
34 """
35 m_copies = [control.copy_data() for control in self.controls]
---> 36 return self.rf.__call__(self.set_local(m_copies, m_array))
37
38 def set_local(self, m, m_array):
~/.local/lib/python3.9/site-packages/pyadjoint/tape.py in wrapper(*args, **kwargs)
107 def wrapper(*args, **kwargs):
108 with stop_annotating():
--> 109 return function(*args, **kwargs)
110
111 return wrapper
~/.local/lib/python3.9/site-packages/pyadjoint/reduced_functional.py in __call__(self, values)
210 range(len(blocks))
211 ):
--> 212 blocks[i].recompute()
213
214 # ReducedFunctional can result in a scalar or an assembled 1-form
~/.local/lib/python3.9/site-packages/pyadjoint/block.py in recompute(self, markings)
345 return
346
--> 347 prepared = self.prepare_recompute_component(inputs, relevant_outputs)
348
349 for idx, out in relevant_outputs:
~/.local/lib/python3.9/site-packages/fenics_adjoint/blocks/solving.py in prepare_recompute_component(self, inputs, relevant_outputs)
440
441 def prepare_recompute_component(self, inputs, relevant_outputs):
--> 442 return self._replace_recompute_form()
443
444 def _replace_recompute_form(self):
~/.local/lib/python3.9/site-packages/fenics_adjoint/blocks/solving.py in _replace_recompute_form(self)
446
447 bcs = self._recover_bcs()
--> 448 lhs = self._replace_form(self.lhs, func=func)
449 rhs = 0
450 if self.linear:
~/.local/lib/python3.9/site-packages/fenics_adjoint/blocks/solving.py in _replace_form(self, form, func)
121 dolfin.Function.assign(func, replace_map[self.func])
122 replace_map[self.func] = func
--> 123 return ufl.replace(form, replace_map)
124
125 def _should_compute_boundary_adjoint(self, relevant_dependencies):
~/.local/lib/python3.9/site-packages/ufl_legacy/algorithms/replace.py in replace(e, mapping)
54 # To fix this would require one to expand derivatives early (which
55 # is not attractive), or make replace lazy too.
---> 56 if has_exact_type(e, CoefficientDerivative):
57 # Hack to avoid circular dependencies
58 from ufl_legacy.algorithms.ad import expand_derivatives
~/.local/lib/python3.9/site-packages/ufl_legacy/algorithms/analysis.py in has_exact_type(a, ufl_type)
86 else:
87 traversal = unique_pre_traversal
---> 88 return any(o._ufl_typecode_ == tc for e in iter_expressions(a) for o in traversal(e))
89
90
~/.local/lib/python3.9/site-packages/ufl_legacy/algorithms/traversal.py in iter_expressions(a)
32 elif isinstance(a, Expr):
33 return (a,)
---> 34 error("Not an UFL type: %s" % str(type(a)))
~/.local/lib/python3.9/site-packages/ufl_legacy/log.py in error(self, *message)
156 "Write error message and raise an exception."
157 self._log.error(*message)
--> 158 raise self._exception_type(self._format_raw(*message))
159
160 def begin(self, *message):
UFLException: Not an UFL type: <class 'ufl.form.Form'>
I understand that the error is in the variational form that I have defined but I have failed do know to get why because the error does not arise when i solve the ode simply. I will very much appreciate any help on this. thank you. The ode equation is \dfrac{dy}{dt} = 2tu(t).