Skip to content

Commit

Permalink
Merge pull request #23 from pariterre/PlayingOnString
Browse files Browse the repository at this point in the history
Playing on string
  • Loading branch information
pariterre committed Oct 30, 2019
2 parents 49ee617 + 888d499 commit 397c7f1
Show file tree
Hide file tree
Showing 14 changed files with 335 additions and 73 deletions.
64 changes: 44 additions & 20 deletions analyses/show_optim_control.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
import numpy as np
import matplotlib.pyplot as plt
import biorbd
from BiorbdViz import BiorbdViz
Expand All @@ -10,9 +11,30 @@
output_files = "AvNPhases"
fun_dyn = utils.dynamics_from_muscles_and_torques
runge_kutta_algo = 'rk45'
nb_nodes = 30
nb_intervals = 30
nb_phases = 2
nb_frame_inter = 500
muscle_plot_mapping = \
[[14, 7, 0, 0, 0, 0], # Trapeze1
[15, 7, 1, 0, 0, 0], # Trapeze2
[16, 8, 0, 0, 0, 0], # Trapeze3
[17, 8, 1, 0, 0, 0], # Trapeze4
[10, 5, 2, 1, 0, 1], # SupraSpin
[8, 5, 0, 1, 0, 1], # InfraSpin
[11, 5, 3, 1, 0, 1], # SubScap
[6, 4, 0, 0, 1, 2], # Pectoral1
[0, 0, 0, 0, 1, 2], # Pectoral2
[1, 0, 1, 0, 1, 2], # Pectoral3
[7, 4, 1, 1, 1, 3], # Deltoid1
[9, 5, 1, 1, 1, 3], # Deltoid2
[2, 1, 0, 1, 1, 3], # Deltoid3
[12, 6, 0, 2, 0, 4], # BicepsLong
[13, 6, 1, 2, 0, 4], # BicepsShort
[3, 2, 0, 2, 1, 5], # TricepsLong
[5, 3, 1, 2, 1, 5], # TricepsMed
[4, 3, 0, 2, 1, 5], # TricepsLat
]
muscle_plot_names = ["Trapèzes", "Coiffe des rotateurs", "Pectoraux", "Deltoïdes", "Biceps", "Triceps"]

# Load the model
m = biorbd.Model(f"../models/{model_name}.bioMod")
Expand All @@ -29,11 +51,11 @@
raise NotImplementedError("Dynamic not implemented yet")

# Read values
t, all_q, all_qdot = utils.read_acado_output_states(f"../optimal_control/Results/States{output_files}.txt", m, nb_nodes,
t, all_q, all_qdot = utils.read_acado_output_states(f"../optimal_control/Results/States{output_files}.txt", m, nb_intervals,
nb_phases)
all_u = utils.read_acado_output_controls(f"../optimal_control/Results/Controls{output_files}.txt", nb_nodes, nb_phases,
all_u = utils.read_acado_output_controls(f"../optimal_control/Results/Controls{output_files}.txt", nb_intervals, nb_phases,
nb_controls)
t_final = utils.organize_time(f"../optimal_control/Results/Parameters{output_files}.txt", t, nb_phases, nb_nodes, parameter=False)
t_final = utils.organize_time(f"../optimal_control/Results/Parameters{output_files}.txt", t, nb_phases, nb_intervals, parameter=False)


# Integrate
Expand Down Expand Up @@ -62,30 +84,32 @@
# plt.plot(t_interp, utils.derive(q_interp, t_interp), '--')
plt.title("Qdot %i" % i)

# for i in range(nb_controls):
# plt.subplot(nb_controls, 3, 3 + (3 * i))
# utils.plot_piecewise_constant(t_final, all_u[i, :])
# plt.title("Acceleration %i" % i)

for i in range(m.nbGeneralizedTorque()):
plt.subplot(m.nbGeneralizedTorque(), 3, 3 + (3 * i))
utils.plot_piecewise_constant(t_final, all_u[m.nbMuscleTotal()+i, :])
plt.title("Torques %i" % i)
plt.tight_layout(w_pad=-1.0, h_pad=-1.0)

# L = []
# for i in range(m.nbMuscleGroups()):
# L.append(m.muscleGroup(i).nbMuscles())
# nb_muscles_max = max(L)
plt.figure("Activations")
cmp = 0
for i in range(m.nbMuscleGroups()):
for j in range(m.muscleGroup(i).nbMuscles()):
#plt.subplot(nb_muscles_max, m.nbMuscleGroups(), i+1+(m.nbMuscleGroups()*j))
plt.subplot(3, 6, cmp+1)
utils.plot_piecewise_constant(t_final, all_u[cmp, :])
plt.title(biorbd.HillType.getRef(m.muscleGroup(i).muscle(j)).name().getString())
if muscle_plot_mapping is None:
for i in range(m.nbMuscleGroups()):
for j in range(m.muscleGroup(i).nbMuscles()):
plt.subplot(3, 6, cmp + 1)
utils.plot_piecewise_constant(t_final, all_u[cmp, :])
plt.title(m.muscleGroup(i).muscle(j).name().getString())
plt.ylim((0, 1))
cmp += 1
else:
nb_row = np.max(muscle_plot_mapping, axis=0)[3] + 1
nb_col = np.max(muscle_plot_mapping, axis=0)[4] + 1
for muscle_map in muscle_plot_mapping:
plt.subplot(nb_row, nb_col, muscle_map[3] * nb_col + muscle_map[4] + 1)
utils.plot_piecewise_constant(t_final, all_u[muscle_map[0], :])
# plt.title(m.muscleGroup(map[1]).muscle(map[2]).name().getString())
plt.title(muscle_plot_names[muscle_map[5]])
plt.ylim((0, 1))
cmp += 1
plt.tight_layout(w_pad=-1.0, h_pad=-1.0)

# plt.ion() # Non blocking plt.show
plt.show()
Expand Down
28 changes: 22 additions & 6 deletions analyses/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,7 @@ def dynamics_from_muscles_and_torques(t_int, states, biorbd_model, u):
for i in range(len(states_dynamics)):
states_dynamics[i] = biorbd.StateDynamics(0, u[i])

biorbd_model.updateMuscles(biorbd_model, states[:nb_q], states[nb_q:], True)
biorbd_model.updateMuscles(states[:nb_q], states[nb_q:], True)
tau = biorbd.Model.muscularJointTorque(biorbd_model, states_dynamics, states[:nb_q], states[nb_q:])

tau_final = tau.get_array() + u[nb_muscle:nb_muscle+nb_tau]
Expand Down Expand Up @@ -178,6 +178,19 @@ def dynamics_from_joint_torque(t_int, states, biorbd_model, u):
return rsh


def dynamics_from_joint_torque_and_contact(t_int, states, biorbd_model, u):
nb_q = biorbd_model.nbQ()
nb_qdot = biorbd_model.nbQdot()
cs = biorbd_model.getConstraints()
qddot = biorbd.Model.ForwardDynamicsConstraintsDirect(biorbd_model, states[:nb_q], states[nb_q:], u, cs).get_array()
rsh = np.ndarray(nb_q + nb_qdot)
for i in range(nb_q):
rsh[i] = states[nb_q+i]
rsh[i + nb_q] = qddot[i]

return rsh


def dynamics_from_accelerations(t_int, states, biorbd_model, u):
nb_q = biorbd_model.nbQ()
nb_qdot = biorbd_model.nbQdot()
Expand Down Expand Up @@ -231,8 +244,13 @@ def integrate_states_from_controls(biorbd_model, t, all_q, all_qdot, all_u, dyn_
q_init = integrated_tp.y[:, -1]
else:
q_init = np.concatenate((all_q[:, interval+1], all_qdot[:, interval+1]))
all_t = np.concatenate((all_t, integrated_tp.t[:-1]))
integrated_state = np.concatenate((integrated_state, integrated_tp.y[:, :-1]), axis=1)

if interval < t.shape[0] - 2:
all_t = np.concatenate((all_t, integrated_tp.t[:-1]))
integrated_state = np.concatenate((integrated_state, integrated_tp.y[:, :-1]), axis=1)
else:
all_t = np.concatenate((all_t, integrated_tp.t[:]))
integrated_state = np.concatenate((integrated_state, integrated_tp.y[:, :]), axis=1)

if verbose:
print(f"Time: {t[interval]}")
Expand Down Expand Up @@ -274,13 +292,11 @@ def plot_piecewise_constant(t, data):
def plot_piecewise_linear(t, data):
plt.plot(t, data)


def derive(q, t):
der = np.ndarray(q.shape)
for i in range(q.shape[1]):
for j in range(q.shape[0]-1):
der[j][i] = (q[j+1][i]-q[j][i])/(t[j+1]-t[j])

return der



19 changes: 11 additions & 8 deletions models/BrasViolon.bioMod
Original file line number Diff line number Diff line change
Expand Up @@ -482,14 +482,17 @@ external_forces 0
position 0 0 0.2425
endmarker

loopconstraint
predecessor Archet
successor Violon
RTpredecessor 0 0 0 xyz 0 0 0.370
RTsuccessor -1.57 0 -2 xyz 0 0 0.2425
axis 1 1 1 1 1 0
stabilizationParameter 0.6
endloopconstraint



// loopconstraint
// predecessor Archet
// successor Violon
// RTpredecessor 0 0 0 xyz 0 0 0.370
// RTsuccessor -1.57 0 -2 xyz 0 0 0.2425
// axis 1 1 1 1 1 0
// stabilizationParameter 0.6
// endloopconstraint

// DEFINITION DES MUSCLES

Expand Down
40 changes: 40 additions & 0 deletions models/ModelTest.bioMod
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
version 3

// Informations générales
root_actuated 1
external_forces 0

// DEFINITION DES SEGMENTS


segment Cube1
translations x //z
// rotations xyz

mass 3.00000

inertia
1.0000000000 0.0000000000 0.0000000000
0.0000000000 1.0000000000 0.0000000000
0.0000000000 0.0000000000 1.0000000000

com 0.0000000000 0.50000000000 0
mesh -1 -1 -1
mesh -1 1 -1
mesh -1 1 1
mesh -1 -1 1
mesh -1 -1 -1
mesh 1 -1 -1
mesh 1 1 -1
mesh -1 1 -1
mesh 1 1 -1
mesh 1 1 1
mesh -1 1 1
mesh 1 1 1
mesh 1 -1 1
mesh -1 -1 1
mesh 1 -1 1
mesh 1 -1 -1
endsegment


1 change: 1 addition & 0 deletions optimal_control/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ target_include_directories(utils PUBLIC
# Write down all files to include here
set(source_files
RepeatedUpAndDownBow.cpp
EocarNphases.cpp
)

# Setup each project
Expand Down
153 changes: 153 additions & 0 deletions optimal_control/EocarNphases.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,153 @@
#include <acado_optimal_control.hpp>
#include <bindings/acado_gnuplot/gnuplot_window.hpp>
#include "BiorbdModel.h"
#include "includes/dynamics.h"
#include "includes/objectives.h"
#include "includes/constraints.h"
#include <vector>

using namespace std;
USING_NAMESPACE_ACADO

/* ---------- Model ---------- */
biorbd::Model m("../../models/ModelTest.bioMod");
unsigned int nQ(m.nbQ()); // states number
unsigned int nQdot(m.nbQdot()); // derived states number
unsigned int nTau(m.nbGeneralizedTorque()); // controls number
unsigned int nMarkers(m.nMarkers()); // markers number
unsigned int nMus(m.nbMuscleTotal()); // muscles number
unsigned int nPhases(2); // phases number
GeneralizedCoordinates Q(nQ), Qdot(nQdot), Qddot(nQdot);
GeneralizedTorque Tau(nTau);
std::vector<std::shared_ptr<biorbd::muscles::StateDynamics>> musclesStates(nMus);

const double t_Start = 0.0;
const double t_End = 4.0;
const int nPoints(31);


int main ()
{
for (unsigned int i=0; i<nMus; ++i)
musclesStates[i] = std::make_shared<biorbd::muscles::StateDynamics>(biorbd::muscles::StateDynamics());
// printf( "nQdot vaut : %d \n", nQdot);
// printf( "nQ vaut : %d \n", nQ);
// printf( "nTau vaut : %d \n", nTau);
// printf( "nMarkers vaut : %d \n", nMarkers);
// printf( "nMus vaut : %d \n", nMus);

/* ---------- INITIALIZATION ---------- */

// DifferentialState x1("", nQ+nQdot, 1);
// Control u1("", nMus+nTau, 1);
// DifferentialEquation f;
// IntermediateState is1("", nQ+nQdot+nMus+nTau, 1);

/* ---------- INITIALIZATION ---------- */
std::vector<DifferentialState> x1;
std::vector<Control> u1;
std::vector<IntermediateState> is1;


for (unsigned int p=0; p<nPhases; ++p){
x1.push_back(DifferentialState("",nQ+nQdot,1));
u1.push_back(Control("", nMus+nTau, 1));
is1.push_back(IntermediateState(nQ + nQdot + nMus + nTau));




for (unsigned int i = 0; i < nQ + nQdot; ++i) // On remplit le IntermediateState avec les X
is1[p](i) = x1[p](i);



for (unsigned int i = 0; i < nMus + nTau; ++i) // Puis avec les controles
is1[p](i+nQ+nQdot) = u1[p](i);

}


// /* ----------- DEFINE OCP ------------- */
// OCP ocp( t_Start , t_End , nPoints-1);
// CFunction mayer(1, mayerVelocity); // Dans une sortie, on fais la somme des vitesses au carré
// CFunction lagrange(1, lagrangeResidualTorques);
// ocp.minimizeLagrangeTerm(lagrange(u1)); // On minimise le terme de Lagrange
// ocp.minimizeMayerTerm(mayer(is1)); // On minimise le terme de Mayer

/* ----------- DEFINE OCP ------------- */
OCP ocp(t_Start, t_End, nPoints);
CFunction lagrangeRT(1, lagrangeResidualTorques);
CFunction lagrangeA(1, lagrangeActivations);
CFunction F(nQ+nQdot, forwardDynamicsFromJointTorque);
DifferentialEquation f ;


/* ------------ CONSTRAINTS ----------- */

for (unsigned int p=0; p<nPhases; ++p){

(f << dot( x1[p] )) == F(is1[p]);


if (p == 0) {

ocp.subjectTo( AT_START, x1[p](0) == 0 );
ocp.subjectTo( AT_END, x1[p](0) == 10 );

ocp.subjectTo( AT_START, x1[p](1) == 0);
ocp.subjectTo( AT_END, x1[p](1) == 0);
}
else {
ocp.subjectTo( 0.0, x1[p], -x1[p-1], 0.0 );
ocp.subjectTo( 0.0, x1[p-1], -x1[p], 0.0 );
}
// ocp.subjectTo( AT_START, x1(2) == 0 );
// ocp.subjectTo( AT_END, x1(2) == 0 );

// ocp.subjectTo( AT_START, x1(3) == 0 );
// ocp.subjectTo( AT_END, x1(3) == 0 );


for (unsigned int i=0; i<nMus; ++i)
ocp.subjectTo(0.01 <= u1[p](i) <= 1);



for (unsigned int i=nMus; i<nMus+nTau; ++i)
ocp.subjectTo(-100 <= u1[p](i) <= 100);
}
ocp.subjectTo( f ) ;
/* ------------ OBJECTIVE ----------- */
Expression sumLagrange = lagrangeRT(u1[0])+ lagrangeA(u1[0]);
for(unsigned int p=1; p<nPhases; ++p)
sumLagrange += lagrangeRT(u1[p]) + lagrangeA(u1[p]);

ocp.minimizeLagrangeTerm( sumLagrange ); // WARNING




/* ---------- VISUALIZATION ------------ */

GnuplotWindow window; // visualize the results in a Gnuplot window
for (unsigned int p=0; p<nPhases; ++p){

window.addSubplot( x1[p](0), "Etat 0 " );
window.addSubplot( x1[p](1), "Etat 1 " );
window.addSubplot( u1[p](0), "CONTROL 1" ) ;
}

/* ---------- OPTIMIZATION ------------ */

OptimizationAlgorithm algorithm( ocp ) ; // construct optimization algorithm ,
algorithm.set(MAX_NUM_ITERATIONS, 500);

algorithm << window;
algorithm.solve();

algorithm.getDifferentialStates("../Results/StatesEocar.txt");
algorithm.getControls("../Results/ControlsEocar.txt");

return 0;
}
Loading

0 comments on commit 397c7f1

Please sign in to comment.