Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def setUp(self):
self.f = ca.MX.sym("f")
self.g = ca.MX.sym("g")
# - m: representing the ball of the mass in kg
# - L: the length of the pendulum bar in meters
# - g: the gravity constant in m/s^2
# - psi: the actuation angle of the manuver in radians, which stays
# constant for this problem
m = 1.0
L = 3.0
g = 9.81
# psi = pl.pi / 2.0
psi = pl.pi / (180.0 * 2)
# System
x = ca.MX.sym("x", 2)
p = ca.MX.sym("p", 1)
u = ca.MX.sym("u", 1)
# f = ca.vertcat([x[1], p[0]/(m*(L**2))*(u-x[0]) - g/L * pl.sin(x[0])])
f = ca.vertcat(x[1], p[0]/(m*(L**2))*(u-x[0]) - g/L * x[0])
phi = x
system = cp.system.System(x = x, u = u, p = p, f = f, phi = phi)
data = pl.loadtxt('data_pendulum.txt')
time_points = data[:500, 0]
numeas = data[:500, 1]
wmeas = data[:500, 2]
N = time_points.size
ydata = pl.array([numeas,wmeas])
udata = [psi] * (N-1)
def __set_cost_function(self, costFunc, mean_ref_s):
""" Define stage cost and terminal cost
"""
# Create GP and cos function symbols
mean_s = ca.MX.sym('mean', self.__Ny)
covar_x_s = ca.MX.sym('covar_x', self.__Ny, self.__Ny)
covar_u_s = ca.MX.sym('covar_u', self.__Nu, self.__Nu)
u_s = ca.MX.sym('u', self.__Nu)
delta_u_s = ca.MX.sym('delta_u', self.__Nu)
Q = ca.MX(self.__Q)
P = ca.MX(self.__P)
R = ca.MX(self.__R)
S = ca.MX(self.__S)
if costFunc is 'quad':
self.__l_func = ca.Function('l', [mean_s, covar_x_s, u_s, covar_u_s, delta_u_s],
[self.__cost_l(mean_s, mean_ref_s, covar_x_s, u_s,
covar_u_s, delta_u_s, Q, R, S)])
self.__lf_func = ca.Function('lf', [mean_s, covar_x_s],
[self.__cost_lf(mean_s, mean_ref_s, covar_x_s, P)])
elif costFunc is 'sat':
self.__l_func = ca.Function('l', [mean_s, covar_x_s, u_s, covar_u_s, delta_u_s],
[self.__cost_saturation_l(mean_s, mean_ref_s,
covar_x_s, u_s, covar_u_s, delta_u_s, Q, R, S)])
self.__lf_func = ca.Function('lf', [mean_s, covar_x_s],
num_hyp = h_ell + h_sf + h_sn + h_m
prior = None
# prior = dict(
# ell_mean = 10,
# ell_std = 10,
# sf_mean = 10,
# sf_std = 10,
# sn_mean = 1e-5,
# sn_std = 1e-5
# )
# Create solver
Y_s = ca.MX.sym('Y', N)
X_s = ca.MX.sym('X', N, Nx)
hyp_s = ca.MX.sym('hyp', 1, num_hyp)
squaredist_s = ca.MX.sym('sqdist', N, N * Nx)
param_s = ca.horzcat(squaredist_s, Y_s)
NLL_func = ca.Function('NLL', [hyp_s, X_s, Y_s, squaredist_s],
[calc_NLL(hyp_s, X_s, Y_s, squaredist_s,
meanFunc=meanFunc, prior=prior)])
nlp = {'x': hyp_s, 'f': NLL_func(hyp_s, X, Y_s, squaredist_s), 'p': param_s}
# NLP solver options
opts = {}
opts['expand'] = True
opts['print_time'] = False
opts['verbose'] = False
opts['ipopt.print_level'] = 1
opts['ipopt.tol'] = 1e-8
opts['ipopt.mu_strategy'] = 'adaptive'
if optimizer_opts is not None:
def construct_upd_z(self, problem=None, lineq_updz=True):
if problem is not None:
self.problem_upd_z = problem
self._lineq_updz = lineq_updz
return 0.
# check if we have linear equality constraints
self._lineq_updz, A, b = self._check_for_lineq()
if not self._lineq_updz:
raise ValueError('For now, only equality constrained QP ' +
'z-updates are allowed!')
x_i = struct_symMX(self.q_i_struct)
x_j = struct_symMX(self.q_ij_struct)
l_i = struct_symMX(self.q_i_struct)
l_ij = struct_symMX(self.q_ij_struct)
t = MX.sym('t')
T = MX.sym('T')
rho = MX.sym('rho')
par = struct_symMX(self.par_struct)
inp = [x_i.cat, l_i.cat, l_ij.cat, x_j.cat, t, T, rho, par.cat]
t0 = t/T
# put symbols in MX structs (necessary for transformation)
x_i = self.q_i_struct(x_i)
x_j = self.q_ij_struct(x_j)
l_i = self.q_i_struct(l_i)
l_ij = self.q_ij_struct(l_ij)
# transform spline variables: only consider future piece of spline
tf = lambda cfs, basis: shift_knot1_fwd(cfs, basis, t0)
self._transform_spline([x_i, l_i], tf, self.q_i)
self._transform_spline([x_j, l_ij], tf, self.q_ij)
# fill in parameters
A = A(par.cat)
self.__Np = Np
self.__clip_negative = clip_negative
""" Create integrator """
# Integrator options
options = {
"abstol" : 1e-6,
"reltol" : 1e-6,
"tf" : dt,
}
if opt is not None:
options.update(opt)
x = ca.MX.sym('x', Nx)
u = ca.MX.sym('u', Nu)
z = ca.MX.sym('z', Nz)
p = ca.MX.sym('p', Np)
par = ca.vertcat(u, p)
dae = {'x': x, 'ode': ode(x,u,z,p), 'p':par}
if alg is not None:
self.__alg0 = ca.Function('alg_0', [x, u],
[alg_0(x, u)])
dae.update({'z':z, 'alg': alg(x, z, u)})
self.Integrator = ca.integrator('Integrator', 'idas', dae, options)
else:
self.Integrator = ca.integrator('Integrator', 'cvodes', dae, options)
#TODO: Fix discrete DAE model
if alg is None:
""" Create discrete RK4 model """
ode_casadi = ca.Function("ode", [x, u, p], [ode(x,u,z,p)])
def _define_mx(self, name, size0, size1, dictionary, value=None):
if value is None:
value = np.zeros((size0, size1))
symbol_name = self._add_label(name)
dictionary[name] = MX.sym(symbol_name, size0, size1)
self._values[name] = value
self.add_to_dict(dictionary[name], name)
return dictionary[name]
mu_n2 = np.zeros((simPoints, number_of_states))
var_n2 = np.zeros((simPoints, number_of_states))
z_n3 = np.concatenate([x0, u])
z_n3.shape = (1, number_of_inputs)
mu3_n = np.zeros((simPoints, number_of_states))
var3_n = np.zeros((simPoints, number_of_states))
covariance2 = np.zeros((number_of_states, number_of_states))
D = number_of_inputs
F = ca.MX.sym('F', npoints, number_of_states)
Xt = ca.MX.sym('X', npoints, number_of_inputs)
hyp = ca.MX.sym('hyp', hyper.shape)
z = ca.MX.sym('z', z_n.shape)
cov = ca.MX.sym('cov', covariance.shape)
cov2 = ca.MX.sym('cov2', 4, 1)
gp_EM = ca.Function('gp', [Xt, F, hyp, z, cov], GP_noisy_input(invK, Xt, F, hyp, D, z, cov))
gp_TA = ca.Function('gp_taylor_approx', [Xt, F, hyp, z, cov2], gp_taylor_approx(invK, Xt, F, hyp, D, z, cov2))
gp_simple = ca.Function('gp_simple', [Xt, F, hyp, z], gp_casadi(invK, hyp, Xt, F, z))
for dt in range(simPoints):
mu, cov = gp_EM.call([X, Y, hyper, z_n, covariance])
mu, cov = mu.full(), cov.full()
mu.shape, cov.shape = (number_of_states), (number_of_states, number_of_states)
mu_n[dt, :], var_n[dt, :] = mu, np.diag(cov)
z_n = ca.vertcat(mu, u).T
covariance[:number_of_states, :number_of_states] = cov
for dt in range(simPoints):
mu, var = gp_simple.call([X, Y, hyper, z_n2])
mu, var = mu.full(), var.full()
def construct_upd_l(self, problem=None):
if problem is not None:
self.problem_upd_l = problem
return 0.
# create parameters
z_ij = struct_symMX(self.q_ij_struct)
l_ij = struct_symMX(self.q_ij_struct)
x_j = struct_symMX(self.q_ij_struct)
t = MX.sym('t')
T = MX.sym('T')
rho = MX.sym('rho')
inp = [x_j, z_ij, l_ij, t, T, rho]
# update lambda
l_ij_new = self.q_ij_struct(l_ij.cat + rho*(x_j.cat - z_ij.cat))
out = [l_ij_new]
# create problem
prob, buildtime = create_function('upd_l_'+str(self._index), inp, out, self.options)
self.problem_upd_l = prob
return buildtime
if len(self.constants) == 0 or not ca.depends_on(equations, constants):
constants = 0
else:
logger.warning('Not all constants have been eliminated. As a result, the affine DAE expression will use a symbolic matrix, as opposed to a numerical sparse matrix.')
if len(self.parameters) == 0 or not ca.depends_on(equations, parameters):
parameters = 0
else:
logger.warning('Not all parameters have been eliminated. As a result, the affine DAE expression will use a symbolic matrix, as opposed to a numerical sparse matrix.')
A = Af(0, constants, parameters)
b = bf(0, constants, parameters)
# Replace veccat'ed states with brand new state vectors so as to avoid the value copy operations induced by veccat.
self._states_vector = ca.MX.sym('states_vector', sum([s.numel() for s in self._symbols(self.states)]))
self._der_states_vector = ca.MX.sym('der_states_vector', sum([s.numel() for s in self._symbols(self.der_states)]))
self._alg_states_vector = ca.MX.sym('alg_states_vector', sum([s.numel() for s in self._symbols(self.alg_states)]))
self._inputs_vector = ca.MX.sym('inputs_vector', sum([s.numel() for s in self._symbols(self.inputs)]))
states_vector = ca.vertcat(self._states_vector, self._der_states_vector, self._alg_states_vector, self._inputs_vector)
equations = [ca.reshape(ca.mtimes(A, states_vector), equations.shape) + b]
setattr(self, equation_list, equations)
if options['expand_mx']:
logger.info("Expanded MX functions will be returned")
self._expand_mx_func = lambda x: x.expand()
logger.info("Finished model simplification")