Skip to content
Snippets Groups Projects
Commit d7c9cf2d authored by Jérémy MACEIRAS's avatar Jérémy MACEIRAS
Browse files

[matlab] Fixed a typo

parent 3e5042fa
No related branches found
No related tags found
No related merge requests found
...@@ -25,20 +25,20 @@ function LQR_infHor ...@@ -25,20 +25,20 @@ function LQR_infHor
param.nbData = 200; %Number of datapoints param.nbData = 200; %Number of datapoints
param.nbRepros = 4; %Number of reproductions param.nbRepros = 4; %Number of reproductions
param.param.nbVarPos = 2; %Dimension of position data (here: x1,x2) param.nbVarPos = 2; %Dimension of position data (here: x1,x2)
param.nbDeriv = 2; %Number of static & dynamic features (D=2 for [x,dx]) param.nbDeriv = 2; %Number of static & dynamic features (D=2 for [x,dx])
param.nbVar = param.param.nbVarPos * param.nbDeriv; %Dimension of state vector in the tangent space param.nbVar = param.nbVarPos * param.nbDeriv; %Dimension of state vector in the tangent space
param.dt = 1E-2; %Time step duration param.dt = 1E-2; %Time step duration
param.rfactor = 4E-2; %Control cost in LQR param.rfactor = 4E-2; %Control cost in LQR
%Control cost matrix %Control cost matrix
R = eye(param.param.nbVarPos) * param.rfactor; R = eye(param.nbVarPos) * param.rfactor;
%Target and desired covariance %Target and desired covariance
param.Mu = [randn(param.param.nbVarPos,1); zeros(param.param.nbVarPos*(param.nbDeriv-1),1)]; param.Mu = [randn(param.nbVarPos,1); zeros(param.nbVarPos*(param.nbDeriv-1),1)];
[Ar,~] = qr(randn(param.param.nbVarPos)); [Ar,~] = qr(randn(param.nbVarPos));
xCov = Ar * diag(rand(param.param.nbVarPos,1)) * Ar' * 1E-1 xCov = Ar * diag(rand(param.nbVarPos,1)) * Ar' * 1E-1
%% Discrete dynamical System settings %% Discrete dynamical System settings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
...@@ -50,25 +50,18 @@ B1d = zeros(param.nbDeriv,1); ...@@ -50,25 +50,18 @@ B1d = zeros(param.nbDeriv,1);
for i=1:param.nbDeriv for i=1:param.nbDeriv
B1d(param.nbDeriv-i+1) = param.dt^i * 1/factorial(i); %Discrete 1D B1d(param.nbDeriv-i+1) = param.dt^i * 1/factorial(i); %Discrete 1D
end end
A = kron(A1d, eye(param.param.nbVarPos)); %Discrete nD A = kron(A1d, eye(param.nbVarPos)); %Discrete nD
B = kron(B1d, eye(param.param.nbVarPos)); %Discrete nD B = kron(B1d, eye(param.nbVarPos)); %Discrete nD
%% discrete LQR with infinite horizon %% discrete LQR with infinite horizon
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q = blkdiag(inv(xCov), zeros(param.param.nbVarPos*(param.nbDeriv-1))); %Precision matrix Q = blkdiag(inv(xCov), zeros(param.nbVarPos*(param.nbDeriv-1))); %Precision matrix
P = solveAlgebraicRiccati_eig_discrete(A, B*(R\B'), (Q+Q')/2); P = solveAlgebraicRiccati_eig_discrete(A, B*(R\B'), (Q+Q')/2);
L = (B' * P * B + R) \ B' * P * A; %Feedback gain (discrete version) L = (B' * P * B + R) \ B' * P * A; %Feedback gain (discrete version)
%Test ratio between kp and kv
if param.nbDeriv>1
kp = eigs(L(:,1:param.param.nbVarPos));
kv = eigs(L(:,param.param.nbVarPos+1:end));
ratio = kv ./ (2 * kp.^.5)
end
for n=1:param.nbRepros for n=1:param.nbRepros
x = [ones(param.param.nbVarPos,1)+randn(param.param.nbVarPos,1)*5E-1; zeros(param.param.nbVarPos*(param.nbDeriv-1),1)]; x = [ones(param.nbVarPos,1)+randn(param.nbVarPos,1)*5E-1; zeros(param.nbVarPos*(param.nbDeriv-1),1)];
for t=1:param.nbData for t=1:param.nbData
r(n).Data(:,t) = x; r(n).Data(:,t) = x;
u = L * (param.Mu - x); %Compute acceleration (with only feedback terms) u = L * (param.Mu - x); %Compute acceleration (with only feedback terms)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment