Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

method to generate an initial guess. Fixes #354 #365

Open
wants to merge 8 commits into
base: master
Choose a base branch
from
Open

method to generate an initial guess. Fixes #354 #365

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
3b87de1
method to generate an initial guess. Fixes #354
Peter230655 Feb 18, 2025
eaa52ed
corrected length of underline
Peter230655 Feb 18, 2025
9d56ba3
corrected as per suggestions
Peter230655 Feb 18, 2025
1a17afb
Add plotting to initial guess test function.
moorepants Feb 19, 2025
f2235d6
corrected error with instance = 0, added two examples
Peter230655 Feb 19, 2025
73e52bb
Merge branch 'master' into linear-initial-guess
Peter230655 Feb 28, 2025
d89338b
removed plot kwarg, added end_time kwarg, expanded the test
Peter230655 Feb 28, 2025
34d2a5d
removed kwarg plot from pendulum
Peter230655 Feb 28, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions examples-gallery/beginner/plot_pendulum_swing_up_fixed_duration.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,9 +78,9 @@
time_symbol=t)

# %%
# Use a random positive initial guess.
initial_guess = np.random.randn(prob.num_free)

# Use a linear guess.
# %%
initial_guess = prob.create_linear_initial_guess(plot=True)
# %%
# Find the optimal solution.
solution, info = prob.solve(initial_guess)
Expand Down
14 changes: 7 additions & 7 deletions examples-gallery/beginner/plot_sliding_block.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
# %%
"""
Block Sliding Over a Hill
=========================
Expand Down Expand Up @@ -174,6 +173,7 @@ def obj_grad(free):
bounds=bounds,
)

initial_guess = prob.create_linear_initial_guess()
solution, info = prob.solve(initial_guess)
solution_list.append(solution)
info_list.append(info)
Expand Down Expand Up @@ -273,15 +273,15 @@ def update(frame):
print(f'Optimal h value is: {solution_list[selection][-1]:.3f}')

# %%
prob_list[selection].plot_objective_value()
_ = prob_list[selection].plot_objective_value()

# %%
# Plot errors in the solution.
prob_list[selection].plot_constraint_violations(solution_list[selection])
_ = prob_list[selection].plot_constraint_violations(solution_list[selection])

# %%
# Plot the trajectories of the block.
prob_list[selection].plot_trajectories(solution_list[selection])
_ = prob_list[selection].plot_trajectories(solution_list[selection])

# %%
# Animate the solution.
Expand All @@ -294,15 +294,15 @@ def update(frame):
print('Message from optimizer:', info_list[selection]['status_msg'])

# %%
prob_list[selection].plot_objective_value()
_ = prob_list[selection].plot_objective_value()

# %%
# Plot errors in the solution.
prob_list[selection].plot_constraint_violations(solution_list[selection])
_ =prob_list[selection].plot_constraint_violations(solution_list[selection])

# %%
# Plot the trajectories of the block.
prob_list[selection].plot_trajectories(solution_list[selection])
_ =prob_list[selection].plot_trajectories(solution_list[selection])

# %%
# Animate the solution.
Expand Down
233 changes: 232 additions & 1 deletion opty/direct_collocation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -244,7 +244,6 @@ def solve(self, free, lagrange=[], zl=[], zu=[]):
gives the status of the algorithm as a message

"""

return super().solve(free, lagrange=lagrange, zl=zl, zu=zu)

def _generate_bound_arrays(self):
Expand Down Expand Up @@ -686,6 +685,238 @@ def parse_free(self, free):

return parse_free(free, n, q, N, variable_duration)

def create_linear_initial_guess(self, plot=False):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not a fan of this plot kwarg. You can very easily do:

guess = problem.create_linear_initial_guess()
problem.plot_trajectories(guess)

It is best not to mix calculation methods with plotting methods.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not a fan of this plot kwarg. You can very easily do:

guess = problem.create_linear_initial_guess()
problem.plot_trajectories(guess)

It is best not to mix calculation methods with plotting methods.

I agree, will correct it.

"""Creates an initial guess that is the linear interpolation between
exact instance constraints. Please see the notes for more information.

Parameters
----------
plot : bool, optional (default=False)
If True, the initial guess will be plotted.

Returns
-------
initial_guess : ndarray shape(n*N + q*N + r + s)
The initial guess for the free variables in the optimization problem.

Notes
-----
- Instance constraints which contain instances of other state variables
are ignored.
- Between successsive instances of a state variable, the values are
interpolated linearly.
- If bounds exist for unknown input trajectories or unknown parameters,
the values are set to the midpoint of the interval of the bounds.
- If bounds exist for the variable time interval, the value is set to
the midpoint of its bounds.
- If one of the range limit of a bound is ``-np.inf`` or ``np.inf``,
the value is set to the other finite limit of bound.
- All else is set to zero.

"""
hilfs = sm.symbols('hilfs')
if self.collocator.instance_constraints is None:
instance_matrix = sm.Matrix([])
else:
instance_matrix = sm.Matrix(self.collocator.instance_constraints)
par_map = self.collocator.known_parameter_map
instance_matrix = instance_matrix.subs({key: par_map[key] for key
in par_map.keys()})
# setting the variable time interval to 1.0 makes getting the node
# numbers of the instance times easy.
instance_matrix = instance_matrix.subs({self.collocator.
node_time_interval: 1.0})

initial_guess = np.zeros(self.num_free)
num_nodes = self.collocator.num_collocation_nodes

# delete instance constraints which contain instances of other state
# variables
row_delete = []
M_new = instance_matrix.copy()
for i, exp1 in enumerate(M_new):
for a in sm.preorder_traversal(exp1):
zaehler = 0
for b in sm.preorder_traversal(a):
if isinstance(b, sm.Function):
zaehler += 1
if zaehler > 1:
row_delete.append(i)
if len(row_delete) > 0:
for row in sorted(row_delete, reverse=True):
M_new.row_del(row)
instance_matrix = M_new

# If the instance constraint is of the from ``x(t) - 0``, or ``x(t)``
# a dummy ``hilfs`` is added so all constraints have the same form.
# This dummy will be set to zero later.
for i, expr in enumerate(instance_matrix):
zaehler = 0
for a in sm.preorder_traversal(expr):
zaehler += 1
if zaehler == 2:
instance_matrix[i, 0] = instance_matrix[i, 0] + hilfs

if self.collocator.instance_constraints is not None:
if not isinstance(self.collocator.node_time_interval, sm.Symbol):
# fixed time interval
liste = []
liste1 = []
for exp1 in instance_matrix:
for a in sm.preorder_traversal(exp1):
liste.append(a)

for entry in liste:
if isinstance(entry, sm.Function):
idx = liste.index(entry)
liste2 = [0, 0, 0, 0]
liste2[0] = entry
liste2[1] = liste[idx + 1]
liste2[2] = liste[idx -1]
liste2[3] = entry.name
liste1.append(liste2)

name_rank = {name.name: i for i, name in enumerate(
self.collocator.state_symbols)}
liste3 = sorted(liste1, key=lambda x: (name_rank[x[3]], x[1]))
duration = self.collocator.node_time_interval * (num_nodes-1)
for state in self.collocator.state_symbols:
state_idx = self.collocator.state_symbols.index(state)
times = []
values = []
for i in liste3:
if state.name == i[3]:
times.append(i[1])
if i[-2] != hilfs:
values.append(-i[2])
else:
values.append(0)
if len(times) == 0:
pass
elif len(times) == 1:
initial_guess[state_idx*num_nodes:(state_idx+1) \
*num_nodes] = values[0]
else:
for i in range(len(times)-1):
start = round(times[i]/duration*num_nodes)
ende = round(times[i+1]/duration*num_nodes)
werte = np.linspace(values[i], values[i+1],
ende-start)
initial_guess[state_idx*num_nodes+start:
state_idx*num_nodes+ende] = werte

# Variable time interval
else:
liste = []
liste1 = []
for exp1 in instance_matrix:
for a in sm.preorder_traversal(exp1):
liste.append(a)
for entry in liste:
if isinstance(entry, sm.Function):
idx = liste.index(entry)
liste2 = [0, 0, 0, 0]
liste2[0] = entry
liste2[1] = liste[idx + 1]
liste2[2] = liste[idx -1]
liste2[3] = entry.name
liste1.append(liste2)

name_rank = {name.name: i for i, name in enumerate(
self.collocator.state_symbols)}
liste3 = sorted(liste1, key=lambda x: (name_rank[x[3]], x[1]))

duration = (num_nodes-1)
for state in self.collocator.state_symbols:
state_idx = self.collocator.state_symbols.index(state)
times = []
values = []
for i in liste3:
if state.name == i[3]:
times.append(i[1])
if i[-2] != hilfs:
values.append(-i[2])
else:
values.append(0)
if len(times) == 0:
pass
elif len(times) == 1:
initial_guess[state_idx*num_nodes:(state_idx+1) \
*num_nodes] = values[0]
else:
for i in range(len(times)-1):
start = int(times[i])
ende = int(times[i+1])
werte = np.linspace(values[i], values[i+1],
ende-start)
initial_guess[state_idx*num_nodes+start:
state_idx*num_nodes+ende] = werte

# set the values of unknown input trajectories.
if self.bounds is not None and \
len(self.collocator.unknown_input_trajectories) > 0:
start_idx = len(self.collocator.state_symbols) * num_nodes
for symb in self.collocator.unknown_input_trajectories:
idx = self.collocator.unknown_input_trajectories.index(symb)
if symb in self.bounds.keys():
if self.bounds[symb][0] == -np.inf:
wert = self.bounds[symb][1]
elif self.bounds[symb][1] == np.inf:
wert = self.bounds[symb][0]
else:
wert = 0.5*(self.bounds[symb][0] + self.bounds[symb][1])
initial_guess[start_idx + idx*num_nodes:start_idx+(idx+1)
*num_nodes] = wert

# set the values of unknown parameters.
if self.bounds is not None and \
len(self.collocator.unknown_parameters) > 0:
start_idx = (len(self.collocator.state_symbols) +
len(self.collocator.unknown_input_trajectories))* num_nodes
for symb in self.collocator.unknown_parameters:
idx = self.collocator.unknown_parameters.index(symb)
if symb in self.bounds.keys():
if self.bounds[symb][0] == -np.inf:
wert = self.bounds[symb][1]
elif self.bounds[symb][1] == np.inf:
wert = self.bounds[symb][0]
else:
wert = 0.5*(self.bounds[symb][0] + self.bounds[symb][1])
initial_guess[start_idx + idx] = wert

# set the value of the variable time interval.
if isinstance(self.collocator.node_time_interval, sm.Symbol):
if self.bounds is not None:
if self.collocator.node_time_interval in self.bounds.keys():
if self.bounds[self.collocator.node_time_interval][0] \
== -np.inf:
wert = self.bounds[self.collocator.
node_time_interval][1]
elif self.bounds[self.collocator.node_time_interval][1] \
== np.inf:
wert = self.bounds[self.collocator.
node_time_interval][0]
else:
wert = 0.5*(self.bounds[self.collocator.
node_time_interval][0] +
self.bounds[self.collocator.
node_time_interval][1])
initial_guess[-1] = wert

if self.bounds is None:
initial_guess[-1] = 1.0 / (num_nodes-1)

elif self.collocator.node_time_interval not in self.bounds.keys():
initial_guess[-1] = 1.0 / (num_nodes-1)
else:
pass

if plot:
self.plot_trajectories(initial_guess)

return initial_guess


class ConstraintCollocator(object):
"""This class is responsible for generating the constraint function and the
sparse Jacobian of the constraint function using direct collocation methods
Expand Down
Loading
Loading