How do I evaluate a symbolic Jacobian at multiple values of the same variable in the same equation?

I am trying to calculate the symbolic Jacobian of a very complicated function, then use the outputs in a separate numerical function by providing inputs with names that match the symbolic variables. A challenge I am running into is that function I am trying to evaluate is built from two different time steps of the same states.

syms L l a r p y Rn Re Li li ai ra pa ya Rei az el azi eli real

where Li, li, ai, ra, pa, ya are second values of the states L, l, a, r, p, y.

%% Transforming from Lla to NED cartesian
E2 = 6.694379990141317e-3; % Eccentricity constant
kappa = 1-E2;
Ce2n = [-sin(L)*cos(l), -sin(L)*sin(l), cos(L); %  ECEF to NED DCM
    -sin(l), cos(l), 0;
    -cos(L)*cos(l), -cos(L)*sin(l), -sin(L)];

p0 = Ce2n*[(Re+a)*cos(L)*cos(l);        %  platform position NED
            (Re+a)*cos(L)*sin(l);
            (kappa*Re+a)*sin(L)];

pj = Ce2n*[(Rei+ai).*cos(Li).*cos(li);  %  emitter positions NED
            (Rei+ai).*cos(Li).*sin(li);
            (kappa*Rei+ai).*sin(Li)];

%% Create the unit LOS vector
cr = cos(r); sr = sin(r);
cp = cos(p); sp = sin(p);
cy = cos(y); sy = sin(y);
        
Cb2n0 = [           cp*cy,             cp*sy,   -sp; ...
        -cr*sy + sr*sp*cy,  cr*cy + sr*sp*sy, sr*cp; ...
        sr*sy + cr*sp*cy,  -sr*cy + cr*sp*sy, cr*cp]';

cr2 = cos(ra); sr2 = sin(ra);
cp2 = cos(pa); sp2 = sin(pa);
cy2 = cos(ya); sy2 = sin(ya);

Cb2nj = [           cp2*cy2,             cp2*sy2,   -sp2; ...
        -cr2*sy2 + sr2*sp2*cy2,  cr2*cy2 + sr2*sp2*sy2, sr2*cp2; ...
        sr2*sy2 + cr2*sp2*cy2,  -sr2*cy2 + cr2*sp2*sy2, cr2*cp2]';

u0 = Cb2n0*[cos(az)*cos(el);
    sin(az)*cos(el);
    sin(el)];
uj = Cb2nj*[cos(azi)*cos(eli);
    sin(azi)*cos(eli);
    sin(eli)];

%% Trilaterate to initialize feature position
p1 = dot(cross(pj-p0,uj),cross(u0,uj))/((norm(cross(u0,uj))^2));
p2 = dot(cross(p0-pj,u0),cross(uj,u0))/((norm(cross(uj,u0))^2)); % pulled 
% from bearings only SLAM for airborne vehicle

feat_pos_n = 1/2*(p0+pj + p1.*u0 + p2.*uj);
feat_pos_g = [L;l;a] + [1/(Rn +a), 0, 0;
                    0, 1/((Re+a)*cos(L)), 0;
                    0, 0, -1]*(feat_pos_n-p0); % pulled from Groves

where feat_pos_g is the final symbolic function I need to take the Jacobian of, and it is messy and untractable to perform the Jacobian by hand. I need the function evaluated at L, l, a, r, p, y and Li, li, ai, ra, pa, ya but since the later aren't different states I can't just do

statevect = [r, p, y, L, l, a, ra, pa, ya, Li, li, ai]

because it would differentiate the sets of state values as separate states, and

statevect = [r, p, y, L, l, a]

won't work because any term with the second set of terms won't be differentiated.

Any help would be greatly appreciated.

[Andrew Relyea]

I believe what I need Jacobian to do is differentiate like terms (r/ra, p/pa/ ...) like they are the same state while keeping track of which value goes with the differentiated term.

[Umar]

Hi @Andrew Relyea,

You mentioned, “I need the function evaluated at L, l, a, r, p, y and Li, li, ai, ra, pa, ya but since the later aren't different states I can't just do

statevect = [r, p, y, L, l, a, ra, pa, ya, Li, li, ai]

because it would differentiate the sets of state values as separate states, and

statevect = [r, p, y, L, l, a]

won't work because any term with the second set of terms won't be differentiated. Any help would be greatly appreciated.”

Please see my response to your comments below.

You are working with a complex system where you have two different time steps represented by two sets of state variables and your goal is to compute the Jacobian of a derived function (feat_pos_g) that combines these states in a specific transformation. However, you are facing a challenge since the symbolic differentiation needs to consider both sets of variables simultaneously. After doing analysis of your code and reviewing the jacobian documentation provided at the link below,

https://www.mathworks.com/help/symbolic/sym.jacobian.html

here is a revised version that correctly constructs the state vector and calculates the Jacobian:

syms L l a r p y Rn Re Li li ai ra pa ya Rei az el azi eli real

% Constants
E2 = 6.694379990141317e-3; % Eccentricity constant
kappa = 1 - E2;

% Transformation from Lla to NED Cartesian
Ce2n = [-sin(L)*cos(l), -sin(L)*sin(l), cos(L); 
      -sin(l), cos(l), 0; 
      -cos(L)*cos(l), -cos(L)*sin(l), -sin(L)];

p0 = Ce2n * [(Re + a) * cos(L) * cos(l); 
            (Re + a) * cos(L) * sin(l); 
            (kappa * Re + a) * sin(L)];

pj = Ce2n * [(Rei + ai) * cos(Li) * cos(li); 
            (Rei + ai) * cos(Li) * sin(li); 
            (kappa * Rei + ai) * sin(Li)];

% Create the unit LOS vector
cr = cos(r); sr = sin(r);
cp = cos(p); sp = sin(p);
cy = cos(y); sy = sin(y);

Cb2n0 = [cp*cy, cp*sy, -sp; 
       -cr*sy + sr*sp*cy, cr*cy + sr*sp*sy, sr*cp; 
       sr*sy + cr*sp*cy, -sr*cy + cr*sp*sy, cr*cp]';

cr2 = cos(ra); sr2 = sin(ra);
cp2 = cos(pa); sp2 = sin(pa);
cy2 = cos(ya); sy2 = sin(ya);

Cb2nj = [cp2*cy2, cp2*sy2, -sp2; 
        -cr2*sy2 + sr2*sp2*cy2, cr2*cy2 + sr2*sp2*sy2, sr2*cp2; 
        sr2*sy2 + cr2*sp2*cy2, -sr2*cy2 + cr2*sp2*sy2, cr2*cp2]';

u0 = Cb2n0 * [cos(az)*cos(el);
             sin(az)*cos(el);
             sin(el)];
uj = Cb2nj * [cos(azi)*cos(eli);
             sin(azi)*cos(eli);
             sin(eli)];

% Trilaterate to initialize feature position
p1 = dot(cross(pj - p0, uj), cross(u0, uj)) / (norm(cross(u0, uj))^2);
p2 = dot(cross(p0 - pj, u0), cross(uj, u0)) / (norm(cross(uj, u0))^2);

feat_pos_n = 1/2 * (p0 + pj + p1 .* u0 + p2 .* uj);

feat_pos_g = [L; l; a] + [1/(Rn + a), 0, 0; 
                        0, 1/((Re + a) * cos(L)), 0; 
                        0, 0, -1] * (feat_pos_n - p0);

V% Create a combined state vector for differentiation
statevect = [L l a r p y Li li ai ra pa ya];

% Calculate the Jacobian of feat_pos_g with respect to the state vector
J_feat_pos_g = jacobian(feat_pos_g, statevect);

% Display the Jacobian
disp(J_feat_pos_g);

Explanation of modified code

Symbolic Variables: define all necessary symbolic variables at once.

Coordinate Transformations: The code maintains your transformations from ECEF to NED Cartesian coordinates.

Unit Vectors: The unit line-of-sight vectors are created using rotation matrices as per your original structure.

Feature Position Calculation: The trilateration method remains unchanged but is now part of a clean computation pipeline.

State Vector Construction: A single combined statevect is constructed from both sets of states ([L l a r p y Li li ai ra pa ya]).This allows MATLAB's jacobian function to differentiate correctly across all required variables.

Jacobian Calculation: Finally, computing the Jacobian using MATLAB's built-in jacobian function and display it.

Please see attached.

Make sure that you validate the outputs from each transformation step for consistency. Now, if performance becomes an issue with more complex systems or larger symbolic expressions, consider numerical approximations or simplifying assumptions where applicable. This approach should meet your requirements and help you effectively navigate through your symbolic differentiation task. If you have any further questions, please let me know.

[Andrew Relyea]

@Umar, thank you for your response. I have reviewed your code and the only differences I can find are cosmetic. I had considered running all twelve variables, but when you do that it creates extra columns in the Jacobian matrix that do no correspond with a state in the original state vector.

I have attached a screen capture to show the form I need the output matrix to take. As you can see, I am using the output to augement another matrix. In the original there is not, for example, a state Li but there is a L. Creating additional columns misaligns the output.

[Umar]

Hi @Andrew Relyea,

You mentioned, “I have attached a screen capture to show the form I need the output matrix to take. As you can see, I am using the output to augement another matrix. In the original there is not, for example, a state Li but there is a L. Creating additional columns misaligns the output.”

Please see my response to your comments below.

After reviewing your attached a screen capture to show the form you need the output matrix to take and reading your comments, please follow my suggestions mentioned below.

Instead of separating your states into two different groups, create a unified state vector that includes both sets of variables. For example:

statevect = [r; p; y; L; l; a; ra; pa; ya; Li; li; ai];

Make sure that your function feat_pos_g is defined in such a way that it takes this unified vector as input. This will allow you to evaluate all states simultaneously:

   function feat_pos_g = calculate_feat_pos_g(statevect)
       % Extract states
       r = statevect(1);
       p = statevect(2);
       y = statevect(3);
       L = statevect(4);
       l = statevect(5);
       a = statevect(6);
       ra = statevect(7);
       pa = statevect(8);
       ya = statevect(9);
       Li = statevect(10);
       li = statevect(11);
       ai = statevect(12);

       % [Insert calculations for feat_pos_n and feat_pos_g here]
   end

Then, use automatic differentiation tools available in MATLAB or numerical differentiation methods to compute the Jacobian. Here’s an example using finite differences:

   function J = calculate_jacobian(statevect)
       epsilon = 1e-6; % Small perturbation for finite difference
       n = length(statevect);
       m = 3; % Assuming feat_pos_g has 3 outputs

       J = zeros(m, n); % Initialize Jacobian

       for i = 1:n
           perturb = zeros(n, 1);
           perturb(i) = epsilon;

           f_plus = calculate_feat_pos_g(statevect + perturb);
           f_minus = calculate_feat_pos_g(statevect - perturb);

           J(:, i) = (f_plus - f_minus) / (2 * epsilon);
       end
   end

When using this Jacobian in subsequent calculations (like augmenting another matrix), you have to make sure that it matches the expected dimensions as indicated in your example matrices. For instance, if you expect a specific structure like Paug, make sure JGp and JGz are appropriately sized.

After implementing these changes, validate your results by comparing against known cases or using test scenarios where the expected outcomes are established.

Hope this helps.

[Umar]

Hi @ Andrew Relyea,

Please let me know if you have made any progress.

[Andrew Relyea]

I have not had a chance to try to implement the numerical option you provided, but I think I see how it could be adapted to meet my needs. I will provide an update after I've run some tests.

[Umar]

Hi @Andrew Relyea,

Thank you for your update regarding the numerical option. I appreciate your willingness to explore its potential adaptation to meet your needs. I look forward to hearing about your findings once you have had the opportunity to conduct your tests. Please do not hesitate to reach out if you have any questions or require further assistance during the implementation process.