Fix TSE and EM propulsion regressions

+PropulsionPkg/Contents.m

@@ -16,4 +16,5 @@
 % PropulsionSizing     - PropulsionPkg.PropulsionSizing is a function.
 % RecomputeSplits      - PropulsionPkg.RecomputeSplits is a function.
 % TestCreatePropArch   - PropulsionPkg.TestCreatePropArch is a function.
+% TestEngineMotorRegressions - PropulsionPkg.TestEngineMotorRegressions is a function.
 % TestPowerAvailable   - PropulsionPkg.TestPowerAvailable is a function.

+PropulsionPkg/PowerSupplementCheck.m

@@ -1,8 +1,8 @@
 function [Psupp] = PowerSupplementCheck(Preq, Arch, Lambda, Eta, TrnType, EtaFan)
 %
-% [Psupp] = PowerSupplementCheck(Preq, Arch, Lambda, Eta, TrnType, EtaFan, itrn)
+% [Psupp] = PowerSupplementCheck(Preq, Arch, Lambda, Eta, TrnType, EtaFan)
 % written by Paul Mokotoff, prmoko@umich.edu
-% last updated: 04 aug 2025
+% last updated: 16 sep 2026
 %
 % In the propulsion architecture, check if any components are either
 % suppling/siphoning power from the gas-turbine engines. If a component is
@@ -100,23 +100,36 @@
         % find the gas-turbine engine being supplemented
         Driving = find((Arch(:, icomp) > 0)' & (TrnType == 1));
 
-        % check if there are multiple (or none) "driving" gas-turbines
-        if (length(Driving) ~= 1)
-
-            % don't know how much each motor is powering each gas-turbine
+        % a shared target needs both an engine and a motor contribution
+        Helping = find((Arch(:, icomp) > 0)' & (TrnType == 0));
+        if (isempty(Driving) || isempty(Helping))
             continue;
-
         end
-
-        % find the electric motors that are supplementing
-        Helping = find((Arch(:, icomp) > 0)' & (TrnType == 0));
-
-        % add the power supplement, accounting for the fan efficiency
-        Psupp(:, Driving) = Psupp(:, Driving) + sum(Preq(:, Helping), 2) .* EtaFan; %#ok<FNDSB>, ignore warning about "find" ... easier to read this way
+
+        % Downstream splits describe fractions of the target demand. Motor
+        % power at the target already includes path losses in Preq, so using
+        % motor input power here would apply the same split a second time.
+        MotorPower = Preq(:, icomp) .* sum(Lambda(icomp, Helping));
+
+        % Only a fan target needs the engine fan-efficiency conversion.
+        TargetEfficiency = 1;
+        if (TrnType(icomp) == 2)
+            TargetEfficiency = EtaFan;
+        end
+
+        % When engines share the target, allocate the motor credit in
+        % proportion to their own contributions to that target.
+        EngineShares = Lambda(icomp, Driving);
+        if (sum(EngineShares) <= 0)
+            EngineShares = ones(size(EngineShares));
+        end
+        Psupp(:, Driving) = Psupp(:, Driving) + ...
+                            MotorPower .* TargetEfficiency .* ...
+                            (EngineShares ./ sum(EngineShares));
 
     end
 end
 
 % ----------------------------------------------------------
 
-end
+end

+PropulsionPkg/ProcessPropArch.m

@@ -2,7 +2,7 @@
 %
 % [Aircraft] = ProcessPropArch(Aircraft)
 % written by Paul Mokotoff, prmoko@umich.edu
-% last updated: 10 mar 2026
+% last updated: 16 sep 2026
 %
 % given a propulsion architecture, find how each gas turbine engine is
 % connected to a propeller. the power required at the propeller is
@@ -55,7 +55,7 @@
 EtaDwn = ones(ncomp, ncomp);
 
 % allocate memory for indexed transmitter arrays
-WhichProp = zeros(1, ntrn);
+WhichProp = repmat({[]}, 1, ntrn);
 HEcoeff   = zeros(1, ntrn);
 
 % create logical array transmitters
@@ -75,53 +75,15 @@
 Preq = zeros(ncomp, 1);
 PSLS = zeros(ncomp, 1);
 
-% loop through all propellers to find indirect gas turbine engines
+% Propagate each propeller's demand through every upstream path. A direct
+% engine connection must not hide a second engine reached through a cable.
 for iprop = PropIdx
-
-    % remember a copy of the index
-    jprop = iprop;
-
-    % remember the index
-    krow = jprop;
-
-    % loop until gas turbine engine is found
-    while (~isempty(krow))
-
-        % check if any are gas turbine engines
-        AnyGTE = logical(sum(Arch(:, krow) > 0 & ID == 1 & itrn, 2));
-
-        if (any(AnyGTE)) % assume only 1 GTE per fan
-
-            % get its index
-            igte = find(AnyGTE);
-
-            % perturb the power vector
-            PowerVector(iprop) = SLSPower(iprop - nsrc);
-
-            % propagate power downstream
-            Pout = PropulsionPkg.PowerFlow(PowerVector, Arch', LamDwn, EtaDwn, -1, 1.0e-06);
-
-            % remember the power output by the engine and the total power
-            % required
-            Preq(igte) = Preq(igte) + Pout(igte);
-            PSLS(igte) = PSLS(igte) + SLSPower(iprop - nsrc);
-
-            % remove the perturbation
-            PowerVector(iprop) = 0;
-
-            % break out of the loop
-            break;
-
-        else
-
-            % remember the indices
-            jrow = krow;
-
-            % search a level deeper
-            [krow, ~] = find(Arch(:, jrow));
-
-        end
-    end    
+    PowerVector(iprop) = SLSPower(iprop - nsrc);
+    Pout = PropulsionPkg.PowerFlow(PowerVector, Arch', LamDwn, EtaDwn, -1, 1.0e-06);
+    ActiveEngines = EngIdx(Pout(EngIdx) > 0);
+    Preq(ActiveEngines) = Preq(ActiveEngines) + Pout(ActiveEngines);
+    PSLS(ActiveEngines) = PSLS(ActiveEngines) + SLSPower(iprop - nsrc);
+    PowerVector(iprop) = 0;
 end
 
 % loop through all propellers to find directly connected gas turbine engines
@@ -134,7 +96,9 @@
     if (~isempty(GTEIdx))
 
        % remember the index
-       WhichProp(GTEIdx - nsrc) = iprop;
+       for igte = GTEIdx(:)'
+           WhichProp{igte - nsrc}(end + 1) = iprop;
+       end
 
     end
 
@@ -150,4 +114,4 @@
 Aircraft.Specs.Propulsion.Engine.HEcoeff = HEcoeff;
 Aircraft.Specs.Propulsion.PropArch.WhichProp = WhichProp;
 
-end
+end

+PropulsionPkg/PropAnalysis.m

@@ -486,13 +486,13 @@
         icol = HasEng(ieng) + nsrc;
 
         % check if it has a propeller
-        iprop = WhichProp(HasEng(ieng));
+        iprop = WhichProp{HasEng(ieng)};
 
         % check if the engine is connected to a propeller
-        if (iprop ~= 0)
+        if (~isempty(iprop))
 
             % get the thrust requirement from the propeller
-            TEng = Tout(ibeg:iend, iprop);
+            TEng = sum(Tout(ibeg:iend, iprop), 2);
 
         else
 
@@ -712,4 +712,4 @@
 
 % ----------------------------------------------------------
 
-end
+end

+PropulsionPkg/PropulsionSizing.m

@@ -195,7 +195,7 @@
             if (Aircraft.Specs.Propulsion.Engine.DesignThrust > 1.0e-06)
 
                 % size the engine
-                Engine = EngineModelPkg.TurbofanNonlinearSizing(Aircraft.Specs.Propulsion.Engine, Psupp(ieng(1)));
+                Engine = EngineModelPkg.TurbofanNonlinearSizing(Aircraft.Specs.Propulsion.Engine, Psupp(ieng(jeng)));
 
                 % turn off engine sizing
                 Engine.Specs.Sizing = 0;                
@@ -207,7 +207,7 @@
                 InletArea = pi * Engine.FanDiam ^ 2 / 4;
 
                 % find the fan connected to the engine
-                ifan = find((Arch(jeng+nsrc, idx) == 1) & Prop);
+                ifan = find((Arch(ieng(jeng)+nsrc, idx) == 1) & Prop);
 
                 % remember the inlet area
                 Aircraft.Specs.Propulsion.InletArea(ifan) = InletArea;
@@ -319,4 +319,4 @@
 
 % ----------------------------------------------------------
 
-end
+end

+PropulsionPkg/RecomputeSplits.m

@@ -2,14 +2,13 @@
 %
 % [Aircraft] = RecomputeSplits(Aircraft, SegBeg, SegEnd)
 % written by Paul Mokotoff, prmoko@umich.edu
-% last updated: 03 sep 2025
+% last updated: 16 sep 2026
 %
 % Re-compute the operational power splits for a "full throttle" setting
 % during the mission.
 %
-% this function previously worked for two elements connected in parallel.
-% it has now been generalized to work for any number of elements connected
-% in parallel.
+% Each parallel target is recomputed from the power actually delivered to
+% that target. Separate targets may share an engine without sharing a split.
 %
 % INPUTS:
 %     Aircraft - structure with information about the aircraft and mission
@@ -53,11 +52,9 @@
 % get the propulsion architecture
 Arch = Aircraft.Specs.Propulsion.PropArch.Arch;
 
-% get the number of parallel connections
-npar = length(ParIndx);
-
 % get the number of sources and transmitters
 nsrc = length(Aircraft.Specs.Propulsion.PropArch.SrcType);
+TrnType = Aircraft.Specs.Propulsion.PropArch.TrnType;
 
 % get the power available (equal to power output for "full throttle" case)
 Pav = Aircraft.Mission.History.SI.Power.Pav(SegBeg:SegEnd, :);
@@ -72,56 +69,64 @@
 % get the number of downstream splits
 nsplit = length(Aircraft.Specs.Power.LamDwn.SLS);
 
-% get a temporary power split
-TmpSplit = LamDwn(1, :);
+OperUps = Aircraft.Specs.Propulsion.PropArch.OperUps;
+OperDwn = Aircraft.Specs.Propulsion.PropArch.OperDwn;
+EtaUps = Aircraft.Specs.Propulsion.PropArch.EtaUps;
+
+% Collect each parallel target's incoming edges. An engine can contribute
+% to several targets, so its helpers cannot be treated as one group.
+Source = [];
+Target = [];
+for itarget = nsrc + (1:length(TrnType))
+    Parents = find(Arch(:, itarget))';
+    Transmitters = Parents(Parents > nsrc & Parents <= nsrc + length(TrnType));
+    Engine = Transmitters(TrnType(Transmitters - nsrc) == 1);
+    Motor = Transmitters(TrnType(Transmitters - nsrc) == 0);
+    if (~isempty(Engine) && ~isempty(Motor))
+        Source = [Source, Parents];
+        Target = [Target, repmat(itarget, 1, length(Parents))];
+    end
+end
 
-% get the original downstream matrix
-OperDwn = PropulsionPkg.EvalSplit(Aircraft.Specs.Propulsion.PropArch.OperDwn, TmpSplit);
+if (isempty(Source) || nsplit == 0)
+    return
+end
 
-% loop through each power split
-for ipar = 1:npar
-
-    % get the index of the main connection
-    imain = ParIndx(ipar);
-
-    % get the supplemental connection(s)
-    isupp = ParConns{imain};
-
-    % account for the source indices
-    imain = imain + nsrc;
-
-    % get all indices
-    jdx = [imain, isupp];
-
-    % find the upstream split
-    iups = find(sum(Arch(jdx, :), 1) == length(isupp) + 1);
-
-    % get the total power output at any given time from those sources
-    Pout = sum(Pav(idx, jdx), 2);
-
-    % find the downstream split contributing
-    for kdx = jdx
-        for isplit = 1:nsplit
-
-            % perturb a split
-            TmpSplit(isplit) = TmpSplit(isplit) + 0.01;
-
-            % get the new matrix
-            OperNew = PropulsionPkg.EvalSplit(Aircraft.Specs.Propulsion.PropArch.OperDwn, TmpSplit);
-
-            % check if the matrices are different
-            if (abs(OperDwn(iups, kdx) - OperNew(iups, kdx)) > 1.0e-06)
-
-                % recompute the power split
-                LamDwn(idx, isplit) = Pav(idx, kdx) ./ Pout;
-
-            end
-
-            % remove the perturbation
-            TmpSplit(isplit) = TmpSplit(isplit) - 0.01;
-
-        end
-    end            
+ncomp = size(Arch, 1);
+DownIndex = sub2ind([ncomp, ncomp], Target, Source);
+UpIndex = sub2ind([ncomp, ncomp], Source, Target);
+for ipoint = find(idx)'
+    Active = Pav(ipoint, Target) > 0;
+    if (~any(Active))
+        continue
+    end
+    Current = LamDwn(ipoint, :);
+    UpMatrix = PropulsionPkg.EvalSplit(OperUps, LamUps(ipoint, :));
+    DownMatrix = PropulsionPkg.EvalSplit(OperDwn, Current);
+    EdgePower = Pav(ipoint, Source(Active))';
+    TargetPower = Pav(ipoint, Target(Active))';
+    UpFraction = UpMatrix(UpIndex(Active));
+    UpEfficiency = EtaUps(UpIndex(Active));
+    Desired = EdgePower .* UpFraction(:) .* UpEfficiency(:) ./ TargetPower;
+    Baseline = DownMatrix(DownIndex(Active));
+    Baseline = Baseline(:);
+    Sensitivity = zeros(sum(Active), nsplit);
+    for isplit = 1:nsplit
+        Perturbed = Current;
+        Perturbed(isplit) = Perturbed(isplit) + 0.01;
+        NewMatrix = PropulsionPkg.EvalSplit(OperDwn, Perturbed);
+        NewFraction = NewMatrix(DownIndex(Active));
+        Sensitivity(:, isplit) = (NewFraction(:) - Baseline) / 0.01;
+    end
+    Change = Sensitivity \ (Desired(:) - Baseline(:));
+    Updated = Current + Change';
+    NewMatrix = PropulsionPkg.EvalSplit(OperDwn, Updated);
+    Actual = NewMatrix(DownIndex(Active));
+    if (any(~isfinite(Updated)) || any(abs(Actual(:) - Desired(:)) > 1.0e-06))
+        error("FAST:EdgePowerInconsistent", ...
+              "Edge power is not consistent. Your power splits are over-constrained.");
+    end
+    LamDwn(ipoint, :) = Updated;
 end
 
 % if any are NaN, return 0 (assume it's from 0 power available)
@@ -132,4 +137,4 @@
 
 % ----------------------------------------------------------
 
-end
+end