Stable outputs when re-running initial conditions #455

rjmccabe3701 · 2021-06-27T04:26:59Z

rjmccabe3701
Jun 27, 2021

I noticed that the FGFDMexec class doesn't reinitialize itself correctly if you call RunIC() multiple times.

Here is a test program that illustrates the weird behavior:

//jsbsim_wrapper.h
#pragma once
struct InitialCondition
{
   double altitude;  //Altitude in meters
   double fpa_gamma; //flight path angle in radians
   double phi_POM;   //Angle (in radians) defining plane of motion

   double phi;       //Euler roll (radians)
   double theta;     //Euler pitch (radians)
   double psi;       //Euler yaw (radians)
   double v_NED[3];  //Velocity in body coords (m/s)
   double pqr[3];    //Euler rates wrt body coords (rad/sec)

   double alpha;     //Angle of attack (radians)

   //Aircraft control surfaces (normalized to [-1,1])
   double elevator;
   double rudder;
   double aileron;

   double throttle;  //Throttle [0,1]
};

struct JsbSimOutputs
{
   //output
   double forces[3];  //Forces (Newtons)
   double moments[3]; //Moments (Newton-meters)
   double loadFactor; //LoadFactor
   double thrust;     //Thrust (Newtons)
   double drag;       //Drag (Newtons)
   double uvw[3];     //Body Vel (m/s)
   double uvw_dot[3]; //Body Accel (m/s^2)
};

void initJsbsim();
void deInitJsbsim();
void calculateOuputs(struct InitialCondition*, struct JsbSimOutputs* outputs);

//jsbsim_wrapper.cpp
#include "jsbsim_wrapper.h"
#include <stdlib.h>
#include <FGFDMExec.h>
#include <FGJSBBase.h>
#include <models/FGPropagate.h>
#include <models/FGAccelerations.h>
#include <models/FGAuxiliary.h>
#include <models/FGMassBalance.h>
#include <models/FGPropulsion.h>
#include <models/FGAerodynamics.h>
#include <models/FGFCS.h>
#include <initialization/FGInitialCondition.h>

namespace {
JSBSim::FGFDMExec* fdmex = nullptr;

const double FeetToMeters = 0.3048;
const double MetersToFeet = 1/FeetToMeters;

const double DegressToRad = 0.0175;
const double RadToDegrees = 1/DegressToRad;

const double PoundsToKg = 0.453592;
const double PoundsToNewtons = 4.44822;

const double dT = 0.05;
}

using std::cout;
using std::endl;

void initJsbsim()
{
   if(fdmex)
      delete fdmex;
   setenv("JSBSIM_DEBUG", "0", 1);
   fdmex = new JSBSim::FGFDMExec;

   fdmex->SetRootDir(SGPath(getenv("JSBSIM_REPO_PATH")));
   fdmex->SetAircraftPath(SGPath("aircraft"));
   fdmex->SetEnginePath(SGPath("aircraft/engine"));

   auto aircraftModel = getenv("JSBSIM_MODEL");
   bool result = fdmex->LoadModel(aircraftModel);
   if(!result)
   {
      cout << "Failed to load model " << aircraftModel << endl;
      return;
   }
}

void deInitJsbsim()
{
   if(fdmex)
      delete fdmex;
   fdmex = nullptr;
   std::cout << __func__ << std::endl;
}

namespace
{
void setInitialConditions(struct InitialCondition* initialCondition)
{
   auto fgic = fdmex->GetIC();

   //No wind
   fgic->SetWindNEDFpsIC(0.0, 0.0, 0.0);

   //Set position
   fgic->SetAltitudeASLFtIC(initialCondition->altitude*MetersToFeet);
   fgic->SetGeodLatitudeDegIC(0);
   fgic->SetLongitudeDegIC(0);

   //Set velocity in NED coords
   fgic->SetVNorthFpsIC(initialCondition->v_NED[0]*MetersToFeet);
   fgic->SetVEastFpsIC(initialCondition->v_NED[1]*MetersToFeet);
   fgic->SetVDownFpsIC(initialCondition->v_NED[2]*MetersToFeet);

   //Set orientation
   fgic->SetPhiDegIC(initialCondition->phi*RadToDegrees);
   fgic->SetThetaDegIC(initialCondition->theta*RadToDegrees);
   fgic->SetPsiDegIC(initialCondition->psi*RadToDegrees);

   //Set flight path angle and Alpha
   fgic->SetFlightPathAngleDegIC(initialCondition->fpa_gamma*RadToDegrees);
   fgic->SetAlphaDegIC(initialCondition->alpha*RadToDegrees);

   fgic->SetPRadpsIC(initialCondition->pqr[0]);
   fgic->SetQRadpsIC(initialCondition->pqr[1]);
   fgic->SetRRadpsIC(initialCondition->pqr[2]);

   //Set control surface commands
   fdmex->GetFCS()->SetDaCmd(initialCondition->aileron);
   fdmex->GetFCS()->SetDeCmd(initialCondition->elevator);
   fdmex->GetFCS()->SetDrCmd(initialCondition->rudder);

   //Throttle
   for (int i = 0; i < fdmex->GetPropulsion()->GetNumEngines(); i++) {
      fdmex->GetFCS()->SetThrottleCmd(i, initialCondition->throttle);
   }
   fdmex->Setdt(dT);

   //Don't care about fuel
   fdmex->GetPropulsion()->SetFuelFreeze(true);

   //All Engines running
   fdmex->GetPropulsion()->InitRunning(-1);
   //Run initial conditions
   if(!fdmex->RunIC())
   {
      cout << "Could not set IC\n";
      return;
   }
}
}

void calculateOuputs(struct InitialCondition* initialCondition,
                                             struct JsbSimOutputs* outputs)
{
   setInitialConditions(initialCondition);
   outputs->loadFactor = fdmex->GetAuxiliary()->GetNlf(); //Assuming units

   auto aero = fdmex->GetAerodynamics();
   outputs->forces[0] = aero->GetForces(1) * PoundsToNewtons;
   outputs->forces[1] = aero->GetForces(2) * PoundsToNewtons;
   outputs->forces[2] = aero->GetForces(3) * PoundsToNewtons;

   outputs->moments[0] = aero->GetMoments(1) * PoundsToNewtons * FeetToMeters;
   outputs->moments[1] = aero->GetMoments(2) * PoundsToNewtons * FeetToMeters;
   outputs->moments[2] = aero->GetMoments(3) * PoundsToNewtons * FeetToMeters;

   outputs->drag = -1* aero->GetForces(1) * PoundsToNewtons;
   outputs->thrust = fdmex->GetPropulsion()->GetForces(1) * PoundsToNewtons;

   //Get Body Velocities
   auto propagate = fdmex->GetPropagate();
   outputs->uvw[0] = propagate->GetUVW(1) * FeetToMeters;
   outputs->uvw[1] = propagate->GetUVW(2) * FeetToMeters;
   outputs->uvw[2] = propagate->GetUVW(3) * FeetToMeters;

   //Get Body Accelerations
   auto accel = fdmex->GetAccelerations();
   outputs->uvw_dot[0] = accel->GetUVWdot(1) * FeetToMeters;
   outputs->uvw_dot[1] = accel->GetUVWdot(2) * FeetToMeters;
   outputs->uvw_dot[2] = accel->GetUVWdot(3) * FeetToMeters;
}

Here is a test program that queries the forces, moments, etc of the aircraft three times with the same
initial conditions.

#include "jsbsim_wrapper.h"
#include <iostream>
#include <string>
#include <stdlib.h>

int main(void)
{
   setenv("JSBSIM_REPO_PATH", "../../../", 1);
   setenv("JSBSIM_MODEL", "f15", 1);
   initJsbsim();
   InitialCondition initCond;
   initCond.altitude = 10000;
   initCond.fpa_gamma = 0;
   initCond.phi_POM = 0;
   initCond.phi = 0;
   initCond.theta = -0.0604;
   initCond.psi = 0;

   initCond.v_NED[0] = 329.4096;
   initCond.v_NED[1] = 0;
   initCond.v_NED[2] = 0;

   initCond.pqr[0] = 0;
   initCond.pqr[1] = 0;
   initCond.pqr[2] = 0;

   initCond.alpha = -0.0604;
   initCond.elevator = 0.2967;
   initCond.rudder = 0;
   initCond.aileron = 0;
   initCond.throttle = 0.8995;

   using std::cout;
   using std::endl;;

   for(int i = 0; i < 3; ++i)
   {
      JsbSimOutputs out;
      calculateOuputs(&initCond, &out);
      cout << "forces = " << out.forces[0] << ", " << out.forces[1] << ", " << out.forces[2]
           << "\nmoments = " << out.moments[0] << ", " << out.moments[1] << ", " <<out.moments[2]
           << "\nloadFactor = " << out.loadFactor << ", thrust = " << out.thrust
           << "\ndrag = " << out.drag
           << "\nuvw = " << out.uvw[0] << ", " << out.uvw[1] << ", " <<out.uvw[2]
           << "\nuvw_dot = " << out.uvw_dot[0] << ", " << out.uvw_dot[1] << ", " <<out.uvw_dot[2]
           << "\n----------------------------\n"
           << endl;
   }
   return 0;
}

I am expecting 3 identical sets of outputs, but it looks like the first time I call it I get different answers:

forces = -42737.5, 0, 256893
moments = 0, 27174.6, 0
loadFactor = -1.80909, thrust = 58810.9
drag = 42737.5
uvw = 328.812, 0, -19.8312
uvw_dot = 1.65332, -1.32186e-17, 26.7752
----------------------------

forces = -42711.8, 0, 257318
moments = 0, 28641.4, 0
loadFactor = -1.72036, thrust = 58810.9
drag = 42711.8
uvw = 328.812, 0, -19.8312
uvw_dot = 1.65502, -1.32186e-17, 26.8034
----------------------------

forces = -42711.8, 0, 257318
moments = 0, 28643, 0
loadFactor = -1.72027, thrust = 58810.9
drag = 42711.8
uvw = 328.812, 0, -19.8312
uvw_dot = 1.65503, -1.32186e-17, 26.8034
----------------------------

What am I doing wrong? Thanks!

seanmcleod · 2021-07-01T09:09:31Z

seanmcleod
Jul 1, 2021
Collaborator

Hmm, so if you look closely enough run 2 and 3 aren't identical either.

One issue I spotted with the F-15 model that may explain these differences is that the control surfaces have <kinematic> elements, so the control position you set via fdmex->GetFCS()->SetDeCmd(initialCondition->elevator) for example won't instantly set the control surface position to the commanded position.

Unless there is specific code during the RunIC() procedure that bypasses the kinematic elements. I haven't double-checked this yet. When a trim function is executed there is code to bypass actuator modeling of lag, maximum rate etc. Actually not sure if kinematic elements are bypassed during trim, need to check.

So as a quick test, remove the kinematic elements from the F-15 model and then see if you get 3 identical results.

2 replies

bcoconni Jul 2, 2021
Maintainer

Actually not sure if kinematic elements are bypassed during trim, need to check.

There is a special logic in FGKinemat.cpp that bypasses the lag introduced by <kinematic> elements:

jsbsim/src/models/flight_control/FGKinemat.cpp

Lines 110 to 113 in a632306

    
           if (fcs->GetTrimStatus()) 
        
             // When trimming the output must be reached in one step 
        
             Output = Input; 
        
           else {

rjmccabe3701 Jul 6, 2021
Author

@seanmcleod Yeah, I noticed this too. I didn't want to make my question more confusing so I left this part out of the example code that I posted:

   if(0)
   {
      //These get translated by the FCS and have an actuation delay.
      //The control surfaces stay 0 if the "flight_control" section is commented out
      //If using the flight_control section, the elevator changes with time
      fdmex->GetFCS()->SetDaCmd(initialCondition->aileron);
      fdmex->GetFCS()->SetDeCmd(initialCondition->elevator);
      fdmex->GetFCS()->SetDrCmd(initialCondition->rudder);
   }
   else
   {
      //Is ignored if the flight_control section exists
      //TODO: these mappings from normalized control surface inputs to degrees comes
      // from the FCS portion of the XML file, but that is commented out to bypass it.
      fdmex->GetFCS()->SetDePos(JSBSim::ofNorm, initialCondition->elevator);
      fdmex->GetFCS()->SetDePos(JSBSim::ofDeg, (initialCondition->elevator*20.5 - 5.5));
      fdmex->GetFCS()->SetDrPos(JSBSim::ofNorm, initialCondition->rudder*30);
      fdmex->GetFCS()->SetDaRPos(JSBSim::ofNorm, initialCondition->aileron*20);
      fdmex->GetFCS()->SetDaLPos(JSBSim::ofNorm, -1*initialCondition->aileron*20);
   }

I am running this with the entire <flight_control name="FCS:f15"> section commented out (it seems that it doesn't let you set the low level controls if FCS is active. Doesn't seem to fix the issue:

forces = -26331.5, 0, 90392.9
moments = 0, 4782.76, 0
loadFactor = -0.653935, thrust = 58810.9
drag = 26331.5
uvw = 328.812, 0, -19.8312
uvw_dot = 2.74177, -1.32186e-17, 15.7289
----------------------------

forces = -26316.4, 0, 90644.1
moments = 0, 5650.27, 0
loadFactor = -0.601456, thrust = 58810.9
drag = 26316.4
uvw = 328.812, 0, -19.8312
uvw_dot = 2.74277, -1.32186e-17, 15.7455
----------------------------

forces = -26316.4, 0, 90644.3
moments = 0, 5651.18, 0
loadFactor = -0.6014, thrust = 58810.9
drag = 26316.4
uvw = 328.812, 0, -19.8312
uvw_dot = 2.74277, -1.32186e-17, 15.7455

seanmcleod · 2021-07-01T14:06:19Z

seanmcleod
Jul 1, 2021
Collaborator

Hmm. on closer inspection of the F15 xml I see the kinematic element is only linked to fcs/elevator-pos-norm and it isn't linked to fcs/elevator-pos-rad which is used in the force and moment functions.

        <channel name="Pitch">
            <summer name="Pitch Trim Sum">
                <input>fcs/elevator-cmd-norm</input>
                <input>fcs/pitch-trim-cmd-norm</input>
                <clipto>
                    <min>-1</min>
                    <max>1</max>
                </clipto>
            </summer>

             <kinematic name="elevator position">
                <input>fcs/pitch-trim-sum</input>
                <traverse>
                    <setting>
                        <position>-1</position>
                        <time>0.6</time>
                    </setting>
                    <setting>
                        <position>1</position>
                        <time>0.6</time>
                    </setting>
                </traverse>
                <output>fcs/elevator-pos-norm</output>
            </kinematic>

            <aerosurface_scale name="elevator position">
                <input>fcs/pitch-trim-sum</input>
                <gain>0.01745</gain>
                <domain>
                    <min>-1</min>
                    <max>1</max>
                </domain>
                <range>
                    <min>-26</min>
                    <max>15</max>
                </range>
                <output>fcs/elevator-pos-rad</output>
            </aerosurface_scale>

And double-checking.

fdm['fcs/elevator-cmd-norm'] = 1

# Initialize the aircraft with initial conditions
fdm.run_ic() 

for i in range(0, 80):
    cmd = fdm['fcs/elevator-cmd-norm']
    norm = fdm['fcs/elevator-pos-norm']
    pos = fdm['fcs/elevator-pos-rad']
    print(cmd, norm, pos)
    fdm.run()

Sure enough fcs/elevator-pos-rad is set immediately, only fcs/elevator-pos-norm lags.

1.0 0.05555555555555556 0.26175
1.0 0.08333333333333334 0.26175
1.0 0.11111111111111112 0.26175
1.0 0.1388888888888889 0.26175
1.0 0.16666666666666669 0.26175
1.0 0.19444444444444448 0.26175
1.0 0.22222222222222227 0.26175
1.0 0.25000000000000006 0.26175
1.0 0.27777777777777785 0.26175
1.0 0.30555555555555564 0.26175
1.0 0.3333333333333334 0.26175
1.0 0.3611111111111112 0.26175
1.0 0.388888888888889 0.26175
1.0 0.4166666666666668 0.26175
1.0 0.4444444444444446 0.26175
1.0 0.4722222222222224 0.26175
1.0 0.5000000000000001 0.26175
1.0 0.5277777777777779 0.26175
1.0 0.5555555555555557 0.26175
1.0 0.5833333333333335 0.26175
1.0 0.6111111111111113 0.26175
1.0 0.6388888888888891 0.26175
1.0 0.6666666666666669 0.26175
1.0 0.6944444444444446 0.26175
1.0 0.7222222222222224 0.26175
1.0 0.7500000000000002 0.26175
1.0 0.777777777777778 0.26175
1.0 0.8055555555555558 0.26175
1.0 0.8333333333333336 0.26175
1.0 0.8611111111111114 0.26175
1.0 0.8888888888888892 0.26175
1.0 0.916666666666667 0.26175
1.0 0.9444444444444448 0.26175
1.0 0.9722222222222225 0.26175
1.0 1.0 0.26175
1.0 1.0 0.26175

3 replies

rjmccabe3701 Jul 6, 2021
Author

Hmm, interesting. What is the point of JSBSim independently tracking elevator-pos-rad and elevator-pos-norm? I would think they would just be a static mapping.

seanmcleod Jul 6, 2021
Collaborator

This is the setup I use in the test pilot training sim I developed. In which the mapping between elevator-pos-rad and elevator-pos-norm is a simple mapping.

      <property value="0"> actuators/pitch-fcs-delay </property>
      <property value="0"> actuators/elevator-lag </property>
      <property value="0.08"> actuators/elevator-rate-limit </property>

        <channel name="Pitch">

            <summer name="Pitch Trim Sum">
              <input>fcs/elevator-cmd-norm</input>
              <input>fcs/pitch-trim-cmd-norm</input>
              <clipto>
                <min>-1</min>
                <max>1</max>
              </clipto>
            </summer>

            <aerosurface_scale name="Elevator Control">
              <input>fcs/pitch-trim-sum</input>
              <range>
                <min>-0.3</min>
                <max>0.3</max>
              </range>
            </aerosurface_scale>

            <dynamic_actuator name="fcs/elevator-actuator">
              <input> fcs/elevator-control </input>
              <delay> actuators/pitch-fcs-delay </delay>
              <lag> actuators/elevator-lag </lag>
              <rate_limit> actuators/elevator-rate-limit </rate_limit>
              <clipto>
                <min> -0.3 </min>
                <max>  0.3 </max>
              </clipto>
              <output>fcs/elevator-pos-rad</output>
            </dynamic_actuator>

            <aerosurface_scale name="Elevator Normalized">
              <input>fcs/elevator-pos-rad</input>
              <domain>
                <min>-0.3</min>
                <max>0.3</max>
              </domain>
              <range>
                <min>-1</min>
                <max>1</max>
              </range>
              <output>fcs/elevator-pos-norm</output>
            </aerosurface_scale>

        </channel>

seanmcleod Jul 6, 2021
Collaborator

@rjmccabe3701 the current F15 FCS is a bit of a 'mess' in terms of how a kinematic element has been applied.

seanmcleod · 2021-07-01T19:14:37Z

seanmcleod
Jul 1, 2021
Collaborator

Implemented a quick python version to see if I could reproduce what you're seeing. Haven't tried matching your exact initial conditions yet.

import jsbsim

# The path supplied to FGFDMExec is the location of the folders "aircraft", "engines" and "systems"
fdm = jsbsim.FGFDMExec('..\\') 

# Load the aircraft 
fdm.load_model('f15') 

# Set engines running
fdm['propulsion/engine[0]/set-running'] = 1
fdm['propulsion/engine[1]/set-running'] = 1

results = []

fdm['ic/h-sl-ft'] = 30000
fdm['ic/vc-kts'] = 285
fdm['ic/gamma-deg'] = 0
fdm['fcs/elevator-cmd-norm'] = -0.05

results = { 'attitude' : [], 'bodyVelocities': [], 'bodyRates': [], 'forces': [], 'moments': [] }

for i in range(0, 10):
    fdm.run_ic()

    attitude = "Attitude: {0}, {1}, {2}".format(fdm['attitude/theta-deg'], fdm['attitude/phi-deg'], fdm['attitude/psi-deg'])
    bodyVelocities = "Velocity (u,v,w): {0}, {1}, {2}".format(fdm['velocities/u-fps'], fdm['velocities/v-fps'], fdm['velocities/w-fps'])
    bodyRates = "Rates (p,q,r): {0}, {1}, {2}".format(fdm['velocities/p-rad_sec'], fdm['velocities/q-rad_sec'], fdm['velocities/r-rad_sec'])
    forces = "Aero Forces (body): (x,y,z): {0}, {1}, {2}".format(fdm['forces/fbx-aero-lbs'], fdm['forces/fby-aero-lbs'], fdm['forces/fbz-aero-lbs'])
    moments = "Aero Moments (body): (x,y,z): {0}, {1}, {2}".format(fdm['moments/l-aero-lbsft'], fdm['moments/m-aero-lbsft'], fdm['moments/n-aero-lbsft'])

    results['attitude'].append(attitude)
    results['bodyVelocities'].append(bodyVelocities)
    results['bodyRates'].append(bodyRates)
    results['forces'].append(forces)
    results['moments'].append(moments)

for key in results:
    for item in results[key]:
        print(item)
    print('')

In terms of results:

Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0

Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0

Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0

Aero Forces (body): (x,y,z): -2989.237635695297, 0.0, -6432.800167182645
Aero Forces (body): (x,y,z): -2989.237635695297, 0.0, -6398.291844587653
Aero Forces (body): (x,y,z): -2989.237635695297, 0.0, -6398.2470975854985
Aero Forces (body): (x,y,z): -2989.237635695297, 0.0, -6398.247039561983
Aero Forces (body): (x,y,z): -2989.237635695297, 0.0, -6398.247039486743
Aero Forces (body): (x,y,z): -2989.237635695297, 0.0, -6398.247039486646
Aero Forces (body): (x,y,z): -2989.237635695297, 0.0, -6398.247039486646
Aero Forces (body): (x,y,z): -2989.237635695297, 0.0, -6398.247039486646
Aero Forces (body): (x,y,z): -2989.237635695297, 0.0, -6398.247039486646
Aero Forces (body): (x,y,z): -2989.237635695297, 0.0, -6398.247039486646

Aero Moments (body): (x,y,z): 0.0, 11334.51069592438, 0.0
Aero Moments (body): (x,y,z): 0.0, 11724.054850442746, 0.0
Aero Moments (body): (x,y,z): 0.0, 11724.559973027981, 0.0
Aero Moments (body): (x,y,z): 0.0, 11724.560628021318, 0.0
Aero Moments (body): (x,y,z): 0.0, 11724.560628870651, 0.0
Aero Moments (body): (x,y,z): 0.0, 11724.560628871754, 0.0
Aero Moments (body): (x,y,z): 0.0, 11724.560628871754, 0.0
Aero Moments (body): (x,y,z): 0.0, 11724.560628871754, 0.0
Aero Moments (body): (x,y,z): 0.0, 11724.560628871754, 0.0
Aero Moments (body): (x,y,z): 0.0, 11724.560628871754, 0.0

So what's interesting and perplexing is that the aero forces and moments are only identical from run 6 onwards.

0 replies

seanmcleod · 2021-07-01T19:19:59Z

seanmcleod
Jul 1, 2021
Collaborator

Hmm, so if I simply switch from the f15 model to the 737 model in my test script above these are the results I get for the 737, basically the aero forces and moments are identical for all 10 runs.

So it appears there is something different about the f15 model.

Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0
Attitude: 1.2722218725854067e-14, 0.0, 0.0

Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0
Velocity (u,v,w): 749.8178322058461, 0.0, 0.0

Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0
Rates (p,q,r): 0.0, 0.0, 0.0

Aero Forces (body): (x,y,z): -11303.996708792423, 0.0, -57760.53469486948
Aero Forces (body): (x,y,z): -11303.996708792423, 0.0, -57760.53469486948
Aero Forces (body): (x,y,z): -11303.996708792423, 0.0, -57760.53469486948
Aero Forces (body): (x,y,z): -11303.996708792423, 0.0, -57760.53469486948
Aero Forces (body): (x,y,z): -11303.996708792423, 0.0, -57760.53469486948
Aero Forces (body): (x,y,z): -11303.996708792423, 0.0, -57760.53469486948
Aero Forces (body): (x,y,z): -11303.996708792423, 0.0, -57760.53469486948
Aero Forces (body): (x,y,z): -11303.996708792423, 0.0, -57760.53469486948
Aero Forces (body): (x,y,z): -11303.996708792423, 0.0, -57760.53469486948
Aero Forces (body): (x,y,z): -11303.996708792423, 0.0, -57760.53469486948

Aero Moments (body): (x,y,z): 0.0, 24710.507804375156, 0.0
Aero Moments (body): (x,y,z): 0.0, 24710.507804375156, 0.0
Aero Moments (body): (x,y,z): 0.0, 24710.507804375156, 0.0
Aero Moments (body): (x,y,z): 0.0, 24710.507804375156, 0.0
Aero Moments (body): (x,y,z): 0.0, 24710.507804375156, 0.0
Aero Moments (body): (x,y,z): 0.0, 24710.507804375156, 0.0
Aero Moments (body): (x,y,z): 0.0, 24710.507804375156, 0.0
Aero Moments (body): (x,y,z): 0.0, 24710.507804375156, 0.0
Aero Moments (body): (x,y,z): 0.0, 24710.507804375156, 0.0
Aero Moments (body): (x,y,z): 0.0, 24710.507804375156, 0.0

0 replies

seanmcleod · 2021-07-01T19:49:00Z

seanmcleod
Jul 1, 2021
Collaborator

So I started looking at and logging each of the lift forces, i.e. aero/coefficient/CLalpha, aero/coefficient/CLDe, aero/coefficient/CLadot. aero/coefficient/CLq to see how they varied across the 10 runs for the f15.

Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, -1978.8144291296521, -0.0, 994.1024816676367
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, -1978.8144291296521, -0.0, 959.5941590726438
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, -1978.8144291296521, -0.0, 959.5494120704899
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, -1978.8144291296521, -0.0, 959.5493540469741
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, -1978.8144291296521, -0.0, 959.5493539717348
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, -1978.8144291296521, -0.0, 959.5493539716374
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, -1978.8144291296521, -0.0, 959.5493539716374
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, -1978.8144291296521, -0.0, 959.5493539716374
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, -1978.8144291296521, -0.0, 959.5493539716374
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, -1978.8144291296521, -0.0, 959.5493539716374

And it looks like I've found the culprit, i.e. CLadot. And sure enough - aero/alphadot-rad_sec varies across the runs.

Alpha-dot: 0.035647635065063944
Alpha-dot: 0.03441019716177031
Alpha-dot: 0.03440859257388077
Alpha-dot: 0.03440859049320886
Alpha-dot: 0.03440859049051084
Alpha-dot: 0.034408590490507346
Alpha-dot: 0.034408590490507346
Alpha-dot: 0.034408590490507346
Alpha-dot: 0.034408590490507346
Alpha-dot: 0.034408590490507346

So we'll need to investigate why alpha-dot varies across initialization runs.

The 737 model doesn't include aero/coefficient/CLadot

0 replies

seanmcleod · 2021-07-01T20:03:27Z

seanmcleod
Jul 1, 2021
Collaborator

So adot is calculated using UVW and UVWdot.

      alpha = atan2(vAeroUVW(eW), vAeroUVW(eU));
      adot = (vAeroUVW(eU)*in.vUVWdot(eW) - vAeroUVW(eW)*in.vUVWdot(eU))/mUW;

And in your original post you listed:

uvw = 328.812, 0, -19.8312
uvw = 328.812, 0, -19.8312
uvw = 328.812, 0, -19.8312

uvw_dot = 1.65332, -1.32186e-17, 26.7752
uvw_dot = 1.65502, -1.32186e-17, 26.8034
uvw_dot = 1.65503, -1.32186e-17, 26.8034

In other words UVW was identical across initialization runs but UVWdot wasn't.

So the question is why the variance in UVWdot across runs.

0 replies

seanmcleod · 2021-07-01T22:25:23Z

seanmcleod
Jul 1, 2021
Collaborator

So FGFDMExec::RunIC() calls Propagate->InitializeDerivatives();.

jsbsim/src/FGFDMExec.cpp

Lines 644 to 654 in 375f5be

    
           bool FGFDMExec::RunIC(void) 
        
           { 
        
             SuspendIntegration(); // saves the integration rate, dt, then sets it to 0.0. 
        
             Initialize(IC.get()); 
        
             Models[eInput]->InitModel(); 
        
             Models[eOutput]->InitModel(); 
        
             Run(); 
        
             Propagate->InitializeDerivatives(); 
        
             ResumeIntegration(); // Restores the integration rate to what it was.

Which uses the current derivative to initialize a past value deque 5 deep, which explains where the magical 6 runs comes from in terms of things stabilizing after 6 runs.

jsbsim/src/models/FGPropagate.cpp

Lines 188 to 196 in 375f5be

    
           // Initialize the past value deques 
        
           void FGPropagate::InitializeDerivatives() 
        
           { 
        
             VState.dqPQRidot.assign(5, in.vPQRidot); 
        
             VState.dqUVWidot.assign(5, in.vUVWidot); 
        
             VState.dqInertialVelocity.assign(5, VState.vInertialVelocity); 
        
             VState.dqQtrndot.assign(5, VState.vQtrndot); 
        
           }

1 reply

bcoconni Jul 2, 2021
Maintainer

Which uses the current derivative to initialize a past value deque 5 deep, which explains where the magical 6 runs comes from in terms of things stabilizing after 6 runs.

The depth of vectors in InitializeDerivatives are only meaningful to the order of the Adams-Bashforth (AB) scheme you are using. By default, the following integration schemes are used:

jsbsim/src/models/FGPropagate.cpp

Lines 93 to 96 in a632306

    
           integrator_rotational_rate = eRectEuler; 
        
           integrator_translational_rate = eAdamsBashforth2; 
        
           integrator_rotational_position = eRectEuler; 
        
           integrator_translational_position = eAdamsBashforth3;

So I don't think there are any relation between the cycle of 6 runs you are describing and the initialization of the double queued vectors in FGPropagate::InitialiazeDerivatives.

bcoconni · 2021-07-02T21:53:07Z

bcoconni
Jul 2, 2021
Maintainer

I think there is a misunderstanding: FGFDMExec::RunIC() is not designed as a complete reset of the simulation. The purpose of RunIC is nothing more than setting a limited number of parameters (mainly position, orientation, velocities and rotational velocities) to the initial conditions defined by the properties which name starts with ic/. Every other parameter is left to its current state, so when RunIC is executed for the second or third time these parameters are not reset to the value they had when RunIC was executed for the first time. This is very likely why you get different values during the 2nd and 3rd call to RunIC.

If you want a complete reinitialization of JSBSim, you should call FGFDMExec::ResetToInitialConditions instead. It triggers an entirely different set of instructions:

jsbsim/src/FGFDMExec.cpp

Lines 690 to 704 in a632306

    
           void FGFDMExec::ResetToInitialConditions(int mode) 
        
           { 
        
             if (Constructing) return; 
        
             if (mode == 1) Output->SetStartNewOutput(); 
        
             InitializeModels(); 
        
             if (Script) 
        
               Script->ResetEvents(); 
        
             else 
        
               Setsim_time(0.0); 
        
             RunIC(); 
        
           }

You can get the same result by setting the property simulation/reset to any value. In the case where you want the outputs to be reinitialized as well, you should set the property (or call the method) to 1.

1 reply

rjmccabe3701 Jul 6, 2021
Author

@bcoconni I tried this:

void setInitialConditions(struct InitialCondition* initialCondition)
{
   cout << "RESETTING\n";
   fdmex->ResetToInitialConditions(1);
   auto fgic = fdmex->GetIC();
   // Rest of my example code

And I still get the discrepancy:

RESETTING
forces = -26459.7, 0, 88268.6
moments = 0, -2554.09, 0
loadFactor = -1.09777, thrust = 58810.9
drag = 26459.7
uvw = 328.812, 0, -19.8312
uvw_dot = 2.73327, -1.32186e-17, 15.5879
----------------------------

RESETTING
forces = -26316.2, 0, 90647.9
moments = 0, 5663.61, 0
loadFactor = -0.600648, thrust = 58810.9
drag = 26316.2
uvw = 328.812, 0, -19.8312
uvw_dot = 2.74279, -1.32186e-17, 15.7458
----------------------------

RESETTING
forces = -26316.2, 0, 90647.9
moments = 0, 5663.61, 0
loadFactor = -0.600648, thrust = 58810.9
drag = 26316.2
uvw = 328.812, 0, -19.8312
uvw_dot = 2.74279, -1.32186e-17, 15.7458
----------------------------

seanmcleod · 2021-07-02T22:07:41Z

seanmcleod
Jul 2, 2021
Collaborator

Here was my thinking, although I haven't had a chance to set a breakpoint and double-check, it was getting late last night and it's getting late again.

So at the time of the very first call to RunIC() the vectors have 5 entries of 0 in them which are used then during the initial Run() call with a single integration. Immediately after this single Run() the InitialiazeDerivatives() is called which now initializes the 5 entries with the result of the 1 integration step during the processing of the first call to RunIC().

Now when the second call to RunIC() comes in during it's single integration step in the call to Run() it is using a vector not with 5 zeros in it, but rather the output from the first call to RunIC().

So basically the first 5 calls to RunIC() aren't calling Run() with the same initial vector to be used for the one integration step. Since the vector isn't re-initialized to all zeros.

Ah, as I was about to post this I just saw your subsequent reply which mentions the use of ResetToInitialConditions(int) which more than likely will be resetting these vectors. I'll test it out.

And it also means we just got lucky with the 737 model which did have the exact same forces and moments for all 10 runs.

1 reply

bcoconni Jul 3, 2021
Maintainer

Now when the second call to RunIC() comes in during it's single integration step in the call to Run() it is using a vector not with 5 zeros in it, but rather the output from the first call to RunIC().

Still, as I mentioned earlier the higher AB order is 3 for the integration of positions which means that it uses 3 elements among the 5 that VState.dqInertialVelocity contains, meaning that the last 2 values of this vector are plainly ignored (and I am not mentioning the other integration schemes that are using even less elements of their respective double queued vectors). How would that trigger a cycle of 6 runs ?

seanmcleod · 2021-07-03T10:23:38Z

seanmcleod
Jul 3, 2021
Collaborator

@rjmccabe3701 so I tested @bcoconni's suggestion of using ResetToInitialConditions(int) instead of RunIC(), which is called as the last step of ResetToInitialConditions(int).

for i in range(0, 10):
    #fdm.run_ic()
    fdm.reset_to_initial_conditions(1)

And sure enough the f15 model now shows identical forces, moments, etc. for all 10 runs.

Aero Forces (body): (x,y,z): -2837.977475144025, 0.0, -8340.35799301889
Aero Forces (body): (x,y,z): -2837.977475144025, 0.0, -8340.35799301889
Aero Forces (body): (x,y,z): -2837.977475144025, 0.0, -8340.35799301889
Aero Forces (body): (x,y,z): -2837.977475144025, 0.0, -8340.35799301889
Aero Forces (body): (x,y,z): -2837.977475144025, 0.0, -8340.35799301889
Aero Forces (body): (x,y,z): -2837.977475144025, 0.0, -8340.35799301889
Aero Forces (body): (x,y,z): -2837.977475144025, 0.0, -8340.35799301889
Aero Forces (body): (x,y,z): -2837.977475144025, 0.0, -8340.35799301889
Aero Forces (body): (x,y,z): -2837.977475144025, 0.0, -8340.35799301889
Aero Forces (body): (x,y,z): -2837.977475144025, 0.0, -8340.35799301889

Aero Moments (body): (x,y,z): 0.0, -13585.026738264562, 0.0
Aero Moments (body): (x,y,z): 0.0, -13585.026738264562, 0.0
Aero Moments (body): (x,y,z): 0.0, -13585.026738264562, 0.0
Aero Moments (body): (x,y,z): 0.0, -13585.026738264562, 0.0
Aero Moments (body): (x,y,z): 0.0, -13585.026738264562, 0.0
Aero Moments (body): (x,y,z): 0.0, -13585.026738264562, 0.0
Aero Moments (body): (x,y,z): 0.0, -13585.026738264562, 0.0
Aero Moments (body): (x,y,z): 0.0, -13585.026738264562, 0.0
Aero Moments (body): (x,y,z): 0.0, -13585.026738264562, 0.0
Aero Moments (body): (x,y,z): 0.0, -13585.026738264562, 0.0

Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, 0.0, -0.0, 922.8458783742286
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, 0.0, -0.0, 922.8458783742286
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, 0.0, -0.0, 922.8458783742286
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, 0.0, -0.0, 922.8458783742286
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, 0.0, -0.0, 922.8458783742286
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, 0.0, -0.0, 922.8458783742286
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, 0.0, -0.0, 922.8458783742286
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, 0.0, -0.0, 922.8458783742286
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, 0.0, -0.0, 922.8458783742286
Lift forces (CLalpha, CLDe, CLq, CLadot): 7417.512114644661, 0.0, -0.0, 922.8458783742286

However it's currently not an ideal replacement for what you're trying to achieve. If you look closely you'll see that CLDe is 0 even though I specifically set a non-zero elevator deflection.

fdm['ic/h-sl-ft'] = 30000
fdm['ic/vc-kts'] = 285
fdm['ic/gamma-deg'] = 0
fdm['fcs/elevator-cmd-norm'] = -0.05

That's because as @bcoconni mentioned ResetToInitialConditions(int) is now resetting all state, then calling RunIC() which resets all the ic/ properties to what you have set them and it then executes a single Run() timestep. So properties you have set like the elevator deflection, engine running state, throttle position etc. will be 0 when the initial Run() timestep is executed.

One option is to expand the argument options to ResetToInitialConditions(int) to allow a flag option to say don't execute RunIC(), I as the caller will do that when I'm ready, e.g.

void FGFDMExec::ResetToInitialConditions(int mode)
{
  if (Constructing) return;

  const int START_NEW_OUTPUT = 0x1;
  const int DONT_EXECUTE_RUN_IC = 0x2;

  if (mode & START_NEW_OUTPUT) Output->SetStartNewOutput();

  InitializeModels();

  if (Script)
    Script->ResetEvents();
  else
    Setsim_time(0.0);

  if (!(mode & DONT_EXECUTE_RUN_IC))
    RunIC();
}

So you would then be able to use code along these lines to achieve what you want.

fdm['ic/h-sl-ft'] = 30000
fdm['ic/vc-kts'] = 285
fdm['ic/gamma-deg'] = 0

for i in range(0, 10):
    fdm.reset_to_initial_conditions(3)

    fdm['propulsion/engine[0]/set-running'] = 1
    fdm['propulsion/engine[1]/set-running'] = 1
    fdm['fcs/elevator-cmd-norm'] = -0.05

    fdm.run_ic()

3 replies

rjmccabe3701 Jul 6, 2021
Author

@seanmcleod This seems like the most reasonable solution. I am playing around with it now. If it works out for me I'll create a merge request with the change.

seanmcleod Jul 6, 2021
Collaborator

@rjmccabe3701 there is already a merge request pending for this - #466

bcoconni Jul 10, 2021
Maintainer

FYI, the PR #466 has been merged to the branch master.

seanmcleod · 2021-07-04T14:22:30Z

seanmcleod
Jul 4, 2021
Collaborator

@rjmccabe3701 I noticed an issue with your drag calculation, not related to the issue with multiple initializations.

   auto aero = fdmex->GetAerodynamics();
   outputs->forces[0] = aero->GetForces(1) * PoundsToNewtons;
   outputs->forces[1] = aero->GetForces(2) * PoundsToNewtons;
   outputs->forces[2] = aero->GetForces(3) * PoundsToNewtons;

   outputs->drag = -1* aero->GetForces(1) * PoundsToNewtons;

The drag issue is that GetForces() returns the forces in the body axis, and not the wind axis, so GetForces(1) isn't the drag force.

JSBSim provides the forces and moments in all 3 axes, i.e. body, stability and wind.

jsbsim/src/models/FGAerodynamics.cpp

Lines 567 to 584 in a632306

    
           PropertyManager->Tie("forces/fbx-aero-lbs",  this, eX, (PMF)&FGAerodynamics::GetForces); 
        
           PropertyManager->Tie("forces/fby-aero-lbs",  this, eY, (PMF)&FGAerodynamics::GetForces); 
        
           PropertyManager->Tie("forces/fbz-aero-lbs",  this, eZ, (PMF)&FGAerodynamics::GetForces); 
        
           PropertyManager->Tie("moments/l-aero-lbsft", this, eL, (PMF)&FGAerodynamics::GetMoments); 
        
           PropertyManager->Tie("moments/m-aero-lbsft", this, eM, (PMF)&FGAerodynamics::GetMoments); 
        
           PropertyManager->Tie("moments/n-aero-lbsft", this, eN, (PMF)&FGAerodynamics::GetMoments); 
        
           PropertyManager->Tie("forces/fwx-aero-lbs",  this, eDrag, (PMF)&FGAerodynamics::GetvFw); 
        
           PropertyManager->Tie("forces/fwy-aero-lbs",  this, eSide, (PMF)&FGAerodynamics::GetvFw); 
        
           PropertyManager->Tie("forces/fwz-aero-lbs",  this, eLift, (PMF)&FGAerodynamics::GetvFw); 
        
           PropertyManager->Tie("forces/fsx-aero-lbs",  this, eX, (PMF)&FGAerodynamics::GetForcesInStabilityAxes); 
        
           PropertyManager->Tie("forces/fsy-aero-lbs",  this, eY, (PMF)&FGAerodynamics::GetForcesInStabilityAxes); 
        
           PropertyManager->Tie("forces/fsz-aero-lbs",  this, eZ, (PMF)&FGAerodynamics::GetForcesInStabilityAxes); 
        
           PropertyManager->Tie("moments/roll-stab-aero-lbsft", this, eRoll, (PMF)&FGAerodynamics::GetMomentsInStabilityAxes); 
        
           PropertyManager->Tie("moments/pitch-stab-aero-lbsft", this, ePitch, (PMF)&FGAerodynamics::GetMomentsInStabilityAxes); 
        
           PropertyManager->Tie("moments/yaw-stab-aero-lbsft", this, eYaw, (PMF)&FGAerodynamics::GetMomentsInStabilityAxes); 
        
           PropertyManager->Tie("moments/roll-wind-aero-lbsft", this, eRoll, (PMF)&FGAerodynamics::GetMomentsInWindAxes); 
        
           PropertyManager->Tie("moments/pitch-wind-aero-lbsft", this, ePitch, (PMF)&FGAerodynamics::GetMomentsInWindAxes); 
        
           PropertyManager->Tie("moments/yaw-wind-aero-lbsft", this, eYaw, (PMF)&FGAerodynamics::GetMomentsInWindAxes);

0 replies

seanmcleod · 2021-07-06T22:02:48Z

seanmcleod
Jul 6, 2021
Collaborator

@rjmccabe3701 out of interest what are you using Jsbsim for?

1 reply

rjmccabe3701 Jul 6, 2021
Author

@seanmcleod we are using it for our physics engine for an autonomous UAV demonstration (that's all I can say here about the project). Right now, I'm using JSBSim to obtain an Energy/Maneuverability (E/M) diagram:

Here is an example:

To get this, I'm wrapping JSBsim with MATLAB to find a trim solution and from that solution I determine the excess energy.

This is why I have interest in re-running the RunICs() many times. It is crazy slow if I need to create a new FGFDMExec object for each function call that MATLAB fmincon evaluates (and there is a fmincon optimization for each point in the E/M grid). Luckily your fix seems to do what I want.

seanmcleod · 2021-07-07T07:25:54Z

seanmcleod
Jul 7, 2021
Collaborator

@rjmccabe3701 always interesting to find out about real-world usage of JSBSim.

Talking of using JSBSim to generate performance related graphs coincidentally there was a discussion on this recently - #428

In terms of using MATLAB for the trim calculations any particular reason you used MATLAB as opposed to the JSBSim's trim routines?

Lastly in terms of MATLAB and JSBSim integration there is a current pull request to add support for JSBSim as a Matlab/Simulink S-Function here - #445

How are you integrating MATLAB and JSBSim?

2 replies

rjmccabe3701 Jul 7, 2021
Author

At some of the extreme points in the EM diagram (for example mach = 0.92, psidot = 15 deg/sec). I get this from JSBSim: Sorry, qdot doesn't appear to be trimmable. Granted I am a newb at aerodynamics and flight dynamics (I am mostly a software and math guy), but a really knowledgeable guy on my team said that JSBsim's trim solver was too simple to solve the extreme cases. So I'm just wrapping JSBSim in MATLAB (and calling it with MATLAB's loadlibrary()/calllib() methods) and using MATLAB's optimization toolbox to minimize a trim metric. It seems to work quite well.

I tried something very similar to #445 and was able to wrap JSBSim in a block:

The idea here was to use MATLAB's native nonlinear system trim-solver/linearization to do the work for me. It sorta worked, but I apparently didn't understand matlab's trim-solver well enough (it seemed more complicated than it needed to be), so I decided to roll my own.

We have also integrated JSBSIm into AFSIM as a custom mover, which may or may not be interesting to others (but AFSIM is ITAR, so I won't be able to share nor talk about the code any further).

bcoconni Jul 10, 2021
Maintainer

At some of the extreme points in the EM diagram (for example mach = 0.92, psidot = 15 deg/sec). I get this from JSBSim: Sorry, qdot doesn't appear to be trimmable. Granted I am a newb at aerodynamics and flight dynamics (I am mostly a software and math guy), but a really knowledgeable guy on my team said that JSBsim's trim solver was too simple to solve the extreme cases.

@rjmccabe3701, I'm not surprised by this feedback. JSBSim default trim algorithm is indeed quite basic but still is surprisingly efficient for such a simple algorithm. However we have a track record of trim failures following very minor changes to the IC (see issue #360 for instance) so it's far from being a robust algorithm.

That's also the reason why there has been some development of alternatives to the legacy trim algorithm but they have all been dropped from JSBSim code due to lack of maintenance.

seanmcleod · 2021-07-07T18:46:22Z

seanmcleod
Jul 7, 2021
Collaborator

In terms of some of the extreme points failing to trim it could be related to the default range of alpha values being used by the trim solver in JSBSim, e.g. take a look at the following comment, default appears to be a maximum of 7.5 degrees, but you can change the range.

#428 (comment)

0 replies

seanmcleod · 2021-07-10T17:07:40Z

seanmcleod
Jul 10, 2021
Collaborator

At some of the extreme points in the EM diagram (for example mach = 0.92, psidot = 15 deg/sec). I get this from JSBSim: Sorry, qdot doesn't appear to be trimmable.

@rjmccabe3701 how were you trimming in JSBSim for your EM graphs? Were you using the pullup trim option and setting ic/targetNlf?

0 replies

seanmcleod · 2021-07-10T17:27:49Z

seanmcleod
Jul 10, 2021
Collaborator

Or were you using the turn trim option and setting a bank angle?

void FGTrim::setupTurn(void){
  double g,phi;
  phi = fgic.GetPhiRadIC();
  if( fabs(phi) > 0.001 && fabs(phi) < 1.56 ) {
    targetNlf = 1 / cos(phi);
    g = fdmex->GetInertial()->GetGravity().Magnitude();
    psidot = g*tan(phi) / fgic.GetUBodyFpsIC();
    cout << targetNlf << ", " << psidot << endl;
  }

}

0 replies

Stable outputs when re-running initial conditions #455

Uh oh!

Uh oh!

rjmccabe3701 Jun 27, 2021

Replies: 16 comments · 14 replies

Uh oh!

seanmcleod Jul 1, 2021 Collaborator

Uh oh!

bcoconni Jul 2, 2021 Maintainer

Uh oh!

rjmccabe3701 Jul 6, 2021 Author

Uh oh!

seanmcleod Jul 1, 2021 Collaborator

Uh oh!

rjmccabe3701 Jul 6, 2021 Author

Uh oh!

seanmcleod Jul 6, 2021 Collaborator

Uh oh!

seanmcleod Jul 6, 2021 Collaborator

Uh oh!

seanmcleod Jul 1, 2021 Collaborator

Uh oh!

seanmcleod Jul 1, 2021 Collaborator

Uh oh!

Uh oh!

seanmcleod Jul 1, 2021 Collaborator

Uh oh!

seanmcleod Jul 1, 2021 Collaborator

Uh oh!

seanmcleod Jul 1, 2021 Collaborator

Uh oh!

Uh oh!

bcoconni Jul 2, 2021 Maintainer

Uh oh!

Uh oh!

bcoconni Jul 2, 2021 Maintainer

Uh oh!

rjmccabe3701 Jul 6, 2021 Author

Uh oh!

Uh oh!

seanmcleod Jul 2, 2021 Collaborator

Uh oh!

bcoconni Jul 3, 2021 Maintainer

Uh oh!

seanmcleod Jul 3, 2021 Collaborator

Uh oh!

rjmccabe3701 Jul 6, 2021 Author

Uh oh!

seanmcleod Jul 6, 2021 Collaborator

Uh oh!

bcoconni Jul 10, 2021 Maintainer

Uh oh!

seanmcleod Jul 4, 2021 Collaborator

Uh oh!

seanmcleod Jul 6, 2021 Collaborator

Uh oh!

Uh oh!

rjmccabe3701 Jul 6, 2021 Author

Uh oh!

seanmcleod Jul 7, 2021 Collaborator

Uh oh!

rjmccabe3701 Jul 7, 2021 Author

Uh oh!

Uh oh!

bcoconni Jul 10, 2021 Maintainer

Uh oh!

rjmccabe3701
Jun 27, 2021

Replies: 16 comments 14 replies

seanmcleod
Jul 1, 2021
Collaborator

bcoconni Jul 2, 2021
Maintainer

rjmccabe3701 Jul 6, 2021
Author

seanmcleod
Jul 1, 2021
Collaborator

rjmccabe3701 Jul 6, 2021
Author

seanmcleod Jul 6, 2021
Collaborator

seanmcleod Jul 6, 2021
Collaborator

seanmcleod
Jul 1, 2021
Collaborator

seanmcleod
Jul 1, 2021
Collaborator

seanmcleod
Jul 1, 2021
Collaborator

seanmcleod
Jul 1, 2021
Collaborator

seanmcleod
Jul 1, 2021
Collaborator

bcoconni Jul 2, 2021
Maintainer

bcoconni
Jul 2, 2021
Maintainer

rjmccabe3701 Jul 6, 2021
Author

seanmcleod
Jul 2, 2021
Collaborator

bcoconni Jul 3, 2021
Maintainer

seanmcleod
Jul 3, 2021
Collaborator

rjmccabe3701 Jul 6, 2021
Author

seanmcleod Jul 6, 2021
Collaborator

bcoconni Jul 10, 2021
Maintainer

seanmcleod
Jul 4, 2021
Collaborator

seanmcleod
Jul 6, 2021
Collaborator

rjmccabe3701 Jul 6, 2021
Author

seanmcleod
Jul 7, 2021
Collaborator

rjmccabe3701 Jul 7, 2021
Author

bcoconni Jul 10, 2021
Maintainer