Program cuases an error

[Peter230655]

Dear Timo,
could you try to run it on your machine?
I could not attach it to the email I sent to you.
Thanks, Peter

"""
Heat Diffusion Process with Inequality
======================================

This is example 10.57 from Betts' book "Practical Methods for Optimal Control
Using Nonlinear Programming", 3rd edition, chapter 10: Test Problems.
It deals with the 'discretization' of a PDE.

**States**

- $y_0, .....y_{N-1}$ : state variables

**Specifieds**

- $u_0, u_{}pi$ : control variables

"""

import numpy as np
import sympy as sm
import sympy.physics.mechanics as me
from opty.direct_collocation import Problem
import matplotlib.pyplot as plt

import time

# %%
# Equations of motion.
t = me.dynamicsymbols._t
T = sm.symbols('T', cls=sm.Function)

N = 20
y = list(me.dynamicsymbols(f'y:{N}'))
uy0 = me.dynamicsymbols('uy0')
u0, upi = me.dynamicsymbols('u0 upi')

faktor1, faktor2 = sm.symbols('faktor1 faktor2')
#Parameters
q1 = 1.e-3
q2 = 1.e-3
a = 0.5
b = 0.2
c = 1.0
delta = np.pi/N


def gdt(k, t):
    x = k * np.pi/N
    return (c * (sm.sin(x) * sm.sin( np.pi*t/5) - a) - b).diff(t)

def faktor(k, t, faktor1, faktor2):
    test = c * (sm.sin(k * np.pi/N) * sm.sin(sm.pi*t/5) - a) - b
    hilfs = 1/(1+sm.exp(-50*test)) * faktor2 + 1/(1+sm.exp(50*test)) * faktor1
    return hilfs

eom = sm.Matrix([
    -y[0].diff(t) + uy0,
    -gdt(1, T(t)) + faktor(1, T(t), faktor1, faktor2) * 1/delta**2 * (y[1] - 2*y[0] + u0),
    *[-gdt(i+1, T(t)) + faktor(i, T(t), faktor1, faktor2) * 1/delta**2 * (y[i+1] - 2*y[i] + y[i-1])
            for i in range(1, N-1)],
    -gdt(N-1, T(t)) + faktor(N-1, T(t), faktor1, faktor2) * 1/delta**2 * (upi - 2*y[N-1] + y[N-2])
])

sm.pprint(eom)
print(len(eom))

# %%
# Solve the optimization problem.
t0, tf = 0.0, 5.0
num_nodes = 501
interval_value = (tf - t0)/(num_nodes - 1)

state_symbols = [uy0] +  y
specified_symbols = (u0, upi)

times = np.linspace(t0, tf, num=num_nodes)

# Specify the objective function and form the gradient.
    # Note, that the values for T(t) are found in times
def obj(free):
    value1 = interval_value * (delta/2 + q1) * sum([free[(N+1)*num_nodes+i]**2
                for i in range(num_nodes)])
    value2 = interval_value * (delta/2 + q2) * sum([free[(N+2)*num_nodes+i]**2
                for i in range(num_nodes)])
    value3 = 0
    for i in range(1, N+1):
        value3 += interval_value * delta * sum([free[i*num_nodes+j]**2
                for j in range(num_nodes)])
    return value1 + value2 + value3

def obj_grad(free):
        grad = np.zeros_like(free)
        grad[(N+1)*num_nodes:(N+2)*num_nodes] = (2 * (delta/2 + q1) *
                interval_value * free[(N+1)*num_nodes:(N+2)*num_nodes])
        grad[(N+2)*num_nodes:(N+3)*num_nodes] = (2 * (delta/2 + q2) *
                interval_value * free[(N+2)*num_nodes:(N+3)*num_nodes])
        for i in range(1, N+1):
            grad[i*num_nodes:(i+1)*num_nodes] = (2 * delta * interval_value *
                free[i*num_nodes:(i+1)*num_nodes])
        return grad

# Specify the symbolic instance constraints, as per the example
instance_constraints = (
        *[y[i].func(t0) for i in range(N)],
    )

#bounds1 = {faktor[i]: (0.0, 1.0) for i in range(5, 14)}
#bound2 = {faktor[i]: (1.0, 100.0) for i in range(5)}
#bounds3 = {faktor[i]: (1.0, 100.0) for i in range(14, N)}
#bounds = {**bounds1, **bound2, **bounds3}
bounds = {faktor1: (0.001, 1.0),
        faktor2: (1.0, 200),
}


# Create the optimization problem and set any options.
prob = Problem(obj,
        obj_grad,
        eom,
        state_symbols,
        num_nodes,
        interval_value,
        instance_constraints=instance_constraints,
        known_trajectory_map={T(t): times},
        bounds=bounds,
    )

# prob.add_option('nlp_scaling_method', 'gradient-based')
prob.add_option('max_iter', 3000 )

# Give some rough estimates for the trajectories.
print('prob num_free',prob.num_free)
initial_guess = np.zeros(prob.num_free)
print('obj', obj(initial_guess))
print('obj_grad',obj_grad(initial_guess).shape)

# Find the optimal solution.
for _ in range(1):
    start = time.time()
    solution, info = prob.solve(initial_guess)
    initial_guess = solution
    print(info['status_msg'])
    print(f'Objective value achieved: {info['obj_val']:.4f}, as per the book ' +
            f'it is {4.68159793*1.e-1}, so the error is: '
            f'{(info['obj_val'] - 4.68159793*1.e-1)/(4.68159793*1.e-1)*100:.3f} % ')
    print(f'Time taken for the simulation: {time.time() - start:.2f} s')

# Plot the optimal state and input trajectories.
prob.plot_trajectories(solution)

# Plot the constraint violations.
prob.plot_constraint_violations(solution)

# Plot the objective function as a function of optimizer iteration.
prob.plot_objective_value()

# %%
x = sm.symbols('x')
t1 = sm.symbols('t1')
g = c * (sm.sin(x) * sm.sin(sm.pi*t1/5) - a) - b
g_lam = sm.lambdify((x, t1), g, cse=True)
Delta = []
for i in range(N):
    delta_h = []
    for j in range(num_nodes):
        delta_h.append(solution[(i+1)*num_nodes + j] - g_lam((i+1)*np.pi/N, times[j]))
    Delta.append(delta_h)

fig, ax = plt.subplots(N ,1, figsize=(10, 2*N))
for i in range(N):
    ax[i].plot(times, Delta[i], label=str(i))
    ax[i].axhline(0, color='black')
    ax[i].legend()
ax[0].set_title('Constraint violation, Must be $\\geq 0.0$')
ax[-1].set_xlabel('time [s]')

[Peter230655]

It bwas simply the "sm.pi" issue.