GEKKO optimization with a sinusoidal load on a mass-spring system

Question

I am working with a dynamic mass-spring-damper system with an external, sinusoidal excitation load. The system can be seen in the picture below. For simplicity I have not included the damper in this post, as it is not very important for this issue.

The dynamic system

3 masses are connected by springs with stiffnesses k1 and k2. The outer masses are connected with a spring to fixed setpoints (setpoint 1 and 2) with spring stiffnesses km. The stiffness km is a fixed stiffness, while k1 and k2 are variables to be optimized by this program. The excitation is a sinusoidal function with an amplitude, frequency and phase. The function is time dependent, since the value of a sinusoidal force varies with each time step.

With a sinusoidal excitation, the system masses will eventually obtain a steady-state sinusoidal response. It will result in response amplitudes x1_amp, x2_amp, x3_amp. It also means that the forces felt by the springs k1 and k2 will oscillate with an amplitude F_spring1,amp and F_spring2,amp at steady-state. I assume that the system is at a steady state response when I perform my optimization and attempt to alter the system response.

Optimization model

My objective is split in two parts:

1. The total motion of the masses should not be too large, as this will result in a large load on the km springs. I therefore wish to limit the motions of the outer masses located at local the coordinates x1 and x3 from the picture.

2. I also wish to keep the spring forces that are connecting the masses relatively low, in order to reduce wear and tear on the springs over time.

The total objective function will therefore be a trade-off between these two goals. The trade-off is represented by weights, where the total weight summarizes to 1.0. Based on these objectives and weights, the stiffnesses k1 and k2 will be optimized. The objective function presented in the picture is an idea I have, and I'm open for suggestions to modifications.

Constraints come from the physical limitations. 1) The system must abide by the equation of motion for a dynamical system.

• I am aware that my EOM is not complete, as the mass and spring terms in it should be matrices for a multiple DOF dynamic system, especially when having distinct stiffnesses. 2) The spring forces cannot be infinitely large 3) The masses cannot crash into one another

Attempt and problems so far

I have attempted a (very) limited demo of this problem. I have 1 main concern, however, which is the time varying excitation load. I am struggling to understand how I can properly implement this in GEKKO in a way that works with how I wish to optimize. One way I can think of is to run a time series simulation, and then somehow check which spring stiffnesses that give the best results. Another is the model predictive control, but I'm not really sure if that will work with the functions I am working with. The code I have so far:

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

# Define system parameters
g = 9.81 #gravity constant

#Setting bounds for spring stiffness and spring forces
k_min = 1e-3
k_max = 1e3
F_springMax = 1e4

num_masses = 3
mass = 1.0
damping = 0.05
wavePeriod = 10
omega = 2.0 * np.pi / wavePeriod    #Wave frequency

#Excitation function parameters
#Phase
wavenumber = omega**2 / g
distance = 1                #1m distance between each mass
phi = wavenumber*distance   #Phase depending on distance from origin and wave frequency
excitationAmplitude = 10

#Setting outer setpoints
setpoint1 = -3              #The origin of this system is in mass nr. 1, with x1 = 0 at time = 0
setpoint2 = 5               #The masses are (for now) only point-masses with no volume. SP2 = 2*distance + 3m

#Initial values
x1Init = 0
x2Init = 1*distance
x3Init = 2*distance

#Weights for obj.func
w1 = 0.5
w2 = 0.5


# Set up GEKKO model
m = GEKKO(remote=False)
m.time = np.linspace(0, 50, 1001)

#listXes = np.array([x1,x2,x3])

k = m.Array(m.Var, num_masses-1, lb=k_min, ub=k_max, value = 100)
x = m.Array(m.Var, num_masses)
x_dot = m.Array(m.Var, num_masses, value=0.0)
springForces = m.Array(m.Var, num_masses-1, lb=0,  ub = F_springMax)

# Define variables and their bounds
x[0] = x1 = m.Var(value = 0, lb = setpoint1)
x[1] = x2 = m.Var(value = x2Init, lb = x1)
x[2] = x3 = m.Var(value = x3Init, lb = x2, ub = setpoint2)

# Define equations
for i in range(num_masses):
    
    F_exc = excitationAmplitude * m.sin(omega * m.time - k*i*distance)

    m.Equation(x_dot[i] == (x[i].dt()))

    #I included damping in this equation of motion, but it should not affect the simulator performance
    #THIS is where my MAIN CONCERN is! The F_exc term containing a time series and the other side of the equation doesnt
    m.Equation(x_dot[i].dt() * mass + damping * x_dot[i] + k[i] * x[i] == F_exc)

    if i < num_masses - 1:     #This equation attempts to connect the masses, but I don't think it works
        m.Connection(var1 = x[i], var2 = x[i+1])

#Set springforces
for i in range(num_masses-1):
    m.Equation(springForces[i] == k[i] * (x[i+1]-x[i]))

# Define objective function
m.Obj(w1 * (((x1-x1Init) - setpoint1)**2 + (setpoint2 - (x3 - x3Init))**2)
    + w2 * ((springForces[0])**2 + (springForces[1])**2))

#NOTE: This objective function is not the final version, but I'm not sure how to write the preferred obj.func. in syntax
#This is mainly due to the need of a maximum value of module responses and spring forces


# Set solver options
m.options.IMODE = 6
m.options.NODES = 3

# Solve optimization problem
m.solve(disp=True)

I get an error when it attempts to read the line:

m.Equation(x_dot[i].dt() * mass + damping * x_dot[i] + k[i] * x[i] == F_exc)

telling me: ValueError: operands could not be broadcast together with shapes (1001,) (2,)

It is understandable that the different dimensions do not run well when put in the same equation, but nevertheless I don't know how to work around it.

Connection of nodes / masses is also something I'm not sure how to do properly in GEKKO. I have an attempt ready at 'm.Connection(...)', but it is not complete. But again, my main concern here is how to optimize when my excitation is a time-varying sinusoidal function, and when I need the system to abide by the equation of motion for a dynamic system.

All help is appreciated! :)

Answer 1

Here are a few modifications that may help:

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

# Define system parameters
g = 9.81 #gravity constant

#Setting bounds for spring stiffness and spring forces
k_min = 1e-3
k_max = 1e3
F_springMax = 1e4

num_masses = 3
mass = 1.0
damping = 0.05
wavePeriod = 10
omega = 2.0 * np.pi / wavePeriod    #Wave frequency

#Excitation function parameters
#Phase
wavenumber = omega**2 / g
distance = 1                #1m distance between each mass
phi = wavenumber*distance   #Phase depending on distance from
                            #  origin and wave frequency
excitationAmplitude = 10

#Setting outer setpoints
setpoint1 = -3              #The origin of this system is in mass
                            #  nr. 1, with x1 = 0 at time = 0
setpoint2 = 5               #The masses are (for now) only 
                            #  point-masses with no volume.
                            #   SP2 = 2*distance + 3m

#Initial values
x1Init = 0
x2Init = 1*distance
x3Init = 2*distance

#Weights for obj.func
w1 = 0.5; w2 = 0.5

# Set up GEKKO model
m = GEKKO(remote=True)
m.time = np.linspace(0, 50, 101)

k = m.Array(m.Var, num_masses-1,
            lb=k_min, ub=k_max, value = 100)
x = m.Array(m.Var, num_masses)
x_dot = m.Array(m.Var, num_masses, value=0.0)
springForces = m.Array(m.Var, num_masses-1,
                       lb=0,  ub = F_springMax)

# Define variables and their bounds
x[0] = x1 = m.Var(value = 0, lb = setpoint1)
x[1] = x2 = m.Var(value = x2Init)
x[2] = x3 = m.Var(value = x3Init,ub = setpoint2)

# avoid collisions
m.Equations([x1<=x2,x2<=x3])

# Define equations
F_exc = [None]*num_masses
t = m.Param(m.time)
for i in range(num_masses):
    F_exc[i] = m.Intermediate(excitationAmplitude 
               * m.sin(omega*t - k[min(i,num_masses-2)]*i*distance))
    m.Equation(x_dot[i] == (x[i].dt()))

    #I included damping in this equation of motion,
    #   but it should not affect the simulator performance
    #THIS is where my MAIN CONCERN is! The F_exc term
    #   containing a time series and the other side
    #   of the equation doesn't
    m.Equation(x_dot[i].dt() * mass 
               + damping * x_dot[i] 
               + k[min(i,num_masses-2)] * x[i] == F_exc[i])

    #if i < num_masses - 1:     #This equation attempts
                               #  to connect the masses,
                               #  but I don't think it works
    #    m.Connection(var1 = x[i], var2 = x[i+1])

#Set springforces
for i in range(num_masses-1):
    m.Equation(springForces[i] == k[i] * (x[i+1]-x[i]))

# Define objective function
m.Minimize(w1 * (((x1-x1Init) - setpoint1)**2
           + (setpoint2 - (x3 - x3Init))**2)
           + w2 * ((springForces[0])**2
           + (springForces[1])**2))

#NOTE: This objective function is not the final version,
#  but I'm not sure how to write the preferred obj.func. 
#  in syntax
#This is mainly due to the need of a maximum value of 
#  module responses and spring forces

# Set solver options
m.options.IMODE = 6
m.options.NODES = 3
m.options.SOLVER = 1

# Solve optimization problem
m.solve(disp=True)

Include m.time as a Parameter instead of as a Numpy array. Remove the x[i]==x[i+1] equations otherwise those will be set to be exactly equal. The inequality constraints are included to avoid collisions.

m.Equations([x1<=x2,x2<=x3])

Setting the lower bound during initialization only sets those bounds based on the initial guess values. See the Gekko documentation and Dynamic Optimization course for additional examples.

Answer 2

I could not add a comment, so I will write my suggestions here:

If F_exc has shape (2,) then I assume it is some kind of function over time. Perhaps it would work if F_exc would be discretised, and so describe forces per time step, and so have shape (1001,).

Working in the time domain can be computationally expensive. Perhaps another suggestion is to determine the eigenfrequencies and -shapes, and convert the force vector to the eigendomain as well using the eigenmatrix. In this way, the amplitudes can easily be determined.