Current and Voltage fluctuate with Extreme Resistor

[Nancy-Salpepi]

Test device
Dell and MacBook Air

Operating System
Win10 and macOS 13.1

Browser
FF and Safari

Problem description
For phetsims/qa#900, when a circuit is built with an extreme resistor and high voltage battery, the current and voltage can both fluctuate.
EDIT: With 'Values' checked, I also see the Real Bulb value fluctuate.

This isn't seen in published.

Steps to reproduce

1. On the lab screen, Check 'Real Bulbs' in the Advanced panel.
2. Build a circuit as shown in the video below, starting with the bulb and moving in a clockwise direction.
3. Add Voltmeter.
4. Close the switch.
5. If you don't see fluctuating values, open the switch, increase the battery voltage and then close switch.

Visuals

voltageCurrentFluctuate.mp4

Troubleshooting information: !!!!! DO NOT EDIT !!!!! Name: ‪Circuit Construction Kit: DC‬ URL: https://phet-dev.colorado.edu/html/circuit-construction-kit-dc/1.3.0-dev.20/phet/circuit-construction-kit-dc_all_phet.html Version: 1.3.0-dev.20 2023-02-15 16:13:54 UTC Features missing: applicationcache, applicationcache, touch Flags: pixelRatioScaling User Agent: Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_7) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/110.0.0.0 Safari/537.36 Language: en-US Window: 1374x712 Pixel Ratio: 1.7999999523162842/1 WebGL: WebGL 1.0 (OpenGL ES 2.0 Chromium) GLSL: WebGL GLSL ES 1.0 (OpenGL ES GLSL ES 1.0 Chromium) Vendor: WebKit (WebKit WebGL) Vertex: attribs: 16 varying: 31 uniform: 1024 Texture: size: 16384 imageUnits: 16 (vertex: 16, combined: 32) Max viewport: 16384x16384 OES_texture_float: true Dependencies JSON: {}

[samreid]

Since the realistic bulb resistance is a function of the current, it seems there is some hysteresis in the system. In time step T0, there is R0 resistance which causes Current I0. The changed current causes the resistance to become R1 in the next time step, which changes the current again. The relevant code for the light bulb resistance is here:

circuit-construction-kit-common/js/model/analysis/LinearTransientAnalysis.ts

Lines 158 to 163 in 28f65ba

	// shift by base so at V=0 the log is 1
	resistorAdapter.resistance = LightBulb.REAL_BULB_COLD_RESISTANCE + coefficient * V / logWithBase( V + base, base );
	circuitElement.resistanceProperty.value = resistorAdapter.resistance;
	}
	} );

Is this to be expected? Do we need to try to circumvent the hysteresis? Or just document it? @arouinfar or @matthew-blackman please advise.

[matthew-blackman]

@samreid, @zepumph and I agree that this issue warrants a collaborative effort, since it may result in some deep-level changes to the model.

[matthew-blackman]

While discussing this issue, the question came up as to whether this could be fixed by implementing SPICE. Tagging the issue here to reopen this discussion if we move ahead with SPICE migration #228

[matthew-blackman]

Subject: [PATCH] Running average bulb patch
---
Index: dot/js/RunningAverage.js
IDEA additional info:
Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP
<+>UTF-8
===================================================================
diff --git a/dot/js/RunningAverage.js b/dot/js/RunningAverage.js
--- a/dot/js/RunningAverage.js	(revision 8e08c4b55fb8bd9ee6901f9e53fa76fcdfd67182)
+++ b/dot/js/RunningAverage.js	(date 1676926499726)
@@ -21,16 +21,7 @@
     this.windowSize = windowSize;
 
     // @private {number[]} - Used circularly.
-    this.samples = new Array( windowSize );
-
-    // @private {number} - We add/subtract samples in a circular array pattern using this index.
-    this.sampleIndex = 0;
-
-    // @private {number} - Total sum of the samples within the window (not yet divided by number of samples)
-    this.total = 0;
-
-    // @private {number} - number of samples received so far
-    this.numSamples = 0;
+    this.samples = [];
 
     this.clear();
   }
@@ -41,13 +32,7 @@
    * @public
    */
   clear() {
-    this.total = 0;
-    this.numSamples = 0;
-
-    // Need to clear all of the samples
-    for ( let i = 0; i < this.windowSize; i++ ) {
-      this.samples[ i ] = 0;
-    }
+    this.samples.length = 0;
   }
 
   /**
@@ -57,7 +42,13 @@
    * @returns {number}
    */
   getRunningAverage() {
-    return this.total / this.numSamples;
+    let total = 0;
+
+    this.samples.forEach( sample => {
+      total += sample;
+    } );
+
+    return total / this.samples.length;
   }
 
   /**
@@ -67,7 +58,7 @@
    * @returns {boolean}
    */
   isSaturated() {
-    return this.numSamples >= this.windowSize;
+    return this.samples.length >= this.windowSize;
   }
 
   /**
@@ -80,23 +71,30 @@
   updateRunningAverage( sample ) {
     assert && assert( typeof sample === 'number' && isFinite( sample ) );
 
-    // Limit at the window size
-    this.numSamples = Math.min( this.windowSize, this.numSamples + 1 );
-
-    // Remove the old sample (will be 0 if there was no sample yet, due to clear())
-    this.total -= this.samples[ this.sampleIndex ];
-    assert && assert( isFinite( this.total ) );
+    this.samples.push( sample );
 
-    // Add in the new sample
-    this.total += sample;
-    assert && assert( isFinite( this.total ) );
-
-    // Overwrite in the array and move to the next index
-    this.samples[ this.sampleIndex ] = sample;
-    this.sampleIndex = ( this.sampleIndex + 1 ) % this.windowSize;
+    while ( this.samples.length >= this.windowSize ) {
+      this.samples.shift();
+    }
 
     return this.getRunningAverage();
   }
+
+  /**
+   * Fills the running average with the given sample value
+   * @public
+   *
+   * @param {number} sample
+   */
+  fillRunningAverage( sample ) {
+    assert && assert( typeof sample === 'number' && isFinite( sample ) );
+
+    this.clear();
+
+    for ( let i = 0; i < this.windowSize; i++ ) {
+      this.updateRunningAverage( sample );
+    }
+  }
 }
 
 dot.register( 'RunningAverage', RunningAverage );
Index: circuit-construction-kit-common/js/model/LightBulb.ts
IDEA additional info:
Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP
<+>UTF-8
===================================================================
diff --git a/circuit-construction-kit-common/js/model/LightBulb.ts b/circuit-construction-kit-common/js/model/LightBulb.ts
--- a/circuit-construction-kit-common/js/model/LightBulb.ts	(revision 75fa275d52e7ada6644f0d4558e0436faaa55441)
+++ b/circuit-construction-kit-common/js/model/LightBulb.ts	(date 1676925986239)
@@ -21,6 +21,7 @@
 import PowerDissipatedProperty from './PowerDissipatedProperty.js';
 import optionize from '../../../phet-core/js/optionize.js';
 import IntentionalAny from '../../../phet-core/js/types/IntentionalAny.js';
+import RunningAverage from '../../../dot/js/RunningAverage.js';
 
 // constants
 
@@ -72,6 +73,8 @@
   public readonly resistanceProperty: NumberProperty;
   private readonly viewTypeProperty: Property<CircuitElementViewType>;
 
+  private readonly resistancesRunningAverage: RunningAverage;
+
   public static createAtPosition( startVertex: Vertex,
                                   endVertex: Vertex,
                                   circuit: Circuit,
@@ -139,6 +142,9 @@
 
     // Fill in the chargePathLength
     this.chargePathLength = this.getPathLength();
+
+    this.resistancesRunningAverage = new RunningAverage( 100 );
+    this.resistancesRunningAverage.fillRunningAverage( LightBulb.REAL_BULB_COLD_RESISTANCE );
   }
 
   // Updates the charge path length when the view changes between lifelike/schematic
@@ -243,6 +249,19 @@
     throw new Error( 'exceeded charge path bounds' );
   }
 
+  public addResistanceValue( newResistanceValue: number ): number {
+    const runningAverage = this.resistancesRunningAverage.updateRunningAverage( newResistanceValue );
+    this.resistanceProperty.value = runningAverage;
+    console.log( 'Num samples: ', this.resistancesRunningAverage.numSamples );
+    console.log( 'Running average: ', this.resistancesRunningAverage );
+    return runningAverage;
+  }
+
+  public clearRunningAverageResistance(): void {
+    this.resistancesRunningAverage.clear();
+    this.resistancesRunningAverage.fillRunningAverage( LightBulb.REAL_BULB_COLD_RESISTANCE );
+  }
+
   public static readonly REAL_BULB_COLD_RESISTANCE = 10;
 }
 
Index: circuit-construction-kit-common/js/model/analysis/LinearTransientAnalysis.ts
IDEA additional info:
Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP
<+>UTF-8
===================================================================
diff --git a/circuit-construction-kit-common/js/model/analysis/LinearTransientAnalysis.ts b/circuit-construction-kit-common/js/model/analysis/LinearTransientAnalysis.ts
--- a/circuit-construction-kit-common/js/model/analysis/LinearTransientAnalysis.ts	(revision 75fa275d52e7ada6644f0d4558e0436faaa55441)
+++ b/circuit-construction-kit-common/js/model/analysis/LinearTransientAnalysis.ts	(date 1676925269384)
@@ -157,8 +157,8 @@
         const coefficient = 3;
 
         // shift by base so at V=0 the log is 1
-        resistorAdapter.resistance = LightBulb.REAL_BULB_COLD_RESISTANCE + coefficient * V / logWithBase( V + base, base );
-        circuitElement.resistanceProperty.value = resistorAdapter.resistance;
+        resistorAdapter.resistance = circuitElement.addResistanceValue( LightBulb.REAL_BULB_COLD_RESISTANCE + coefficient * V / logWithBase( V + base, base ) );
+        // circuitElement.resistanceProperty.value = resistorAdapter.resistance;
       }
     } );
 
@@ -197,6 +197,7 @@
       // Clear disconnected isReal light bulbs
       if ( circuitElement instanceof LightBulb && circuitElement.isReal ) {
         circuitElement.resistanceProperty.value = LightBulb.REAL_BULB_COLD_RESISTANCE;
+        circuitElement.clearRunningAverageResistance();
       }
     } );
 

[samreid]

@matthew-blackman and @zepumph and I investigated this. We primarily tried adding a window for smoothing the resistance values of the real bulb. But we were still able to create situations where there was a fluctuation in the readouts. The smoothing window didn't seem to move the fluctuation as far out as we expected, so we aren't sure if we implemented it correctly or if there is a flaw in that logic. We would like to discuss ideas with @arouinfar for how to proceed.

@matthew-blackman also recommends one more idea: Computing the bulb resistance once based on the topology of the circuit and what it is connected to, but only updating it once.

UPDATE: We also saw a former algorithm that used to run one solution with the cold value to determine the operating value. Then runs the circuit solution again with the operating value. We reverted that in 040a20e, possibly (though it's unclear that was working as we expect). More details about non-ohmic behavior and oscillation in #10

[matthew-blackman]

I explored the following code this morning in LinearTransitAnalysis:

const dI = circuitResult.getFinalState().ltaSolution!.getCurrent( resistorAdapter );
const newResistance = 10 + 10 * Math.abs( dI );
circuitElement.resistanceProperty.value += 0.001 * ( newResistance - resistorAdapter.resistance );
resistorAdapter.resistance = circuitElement.resistanceProperty.value;

Setting the resistance of the real bulb based on the current through it doesn't appear to affect this issue, but adding the "drag" to the resistance change with the line circuitElement.resistanceProperty.value += 0.001 * ( newResistance - resistorAdapter.resistance ); appears to have the same affect as the RunningAverage, but with less code.

A major piece of evidence is the fact that bumping the 'Wire Resistivity' up even slightly from the minimum value has a huge effect (factor of 10^6 or more) on the current fluctuation. Could this be a deeper issue within the LTA?

[matthew-blackman]

I believe I have found a viable solution that does not impact any other behavior in either screen.

Setting the minimum wire resistivity to 1E-7 (was originally 1E-10) causes the fluctuations in the real bulb to be undetectable by any sensors or values displays. The rest of the circuit behavior is unchanged in my testing (eg the voltage across all wires is still 0.000 V when the resistivity is set to 'tiny').

This appears to be a workable and low-cost solution. @samreid and @zepumph let's discuss if this is viable for the current RC. We could implement a more robust solution in the future with SPICE, but if this works then I think it's a reasonable stopgap fix for the PhET-iO release.

[matthew-blackman]

We would like to discuss with @arouinfar before moving ahead with possible options.

[samreid]

Patch for discussion:

Subject: [PATCH] Factor out includeExtremeElements, see https://github.com/phetsims/circuit-construction-kit-common/issues/900
---
Index: js/model/analysis/LinearTransientAnalysis.ts
IDEA additional info:
Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP
<+>UTF-8
===================================================================
diff --git a/js/model/analysis/LinearTransientAnalysis.ts b/js/model/analysis/LinearTransientAnalysis.ts
--- a/js/model/analysis/LinearTransientAnalysis.ts	(revision 75fa275d52e7ada6644f0d4558e0436faaa55441)
+++ b/js/model/analysis/LinearTransientAnalysis.ts	(date 1677001486835)
@@ -55,6 +55,8 @@
     const voltageSourceMap = new Map<LTAResistiveBattery, VoltageSource>();
     const capacitorMap = new Map<LTACapacitor, Capacitor>();
     const inductorMap = new Map<LTAInductor, Inductor>();
+
+    let hasRealBulbs = false;
     for ( let i = 0; i < circuit.circuitElements.length; i++ ) {
       const circuitElement = circuit.circuitElements[ i ];
 
@@ -93,6 +95,11 @@
           );
           resistorMap.set( resistorAdapter, circuitElement );
           ltaResistors.push( resistorAdapter );
+
+          if ( circuitElement instanceof LightBulb && circuitElement.isReal ) {
+            resistorAdapter.resistance = LightBulb.REAL_BULB_COLD_RESISTANCE;
+            hasRealBulbs = true;
+          }
         }
         else if ( circuitElement instanceof Switch && !circuitElement.isClosedProperty.value ) {
 
@@ -134,34 +141,40 @@
     }
 
     // Solve the system
-    const ltaCircuit = new LTACircuit( ltaResistors, ltaBatteries, ltaCapacitors, ltaInductors );
-    const circuitResult = ltaCircuit.solveWithSubdivisions( TIMESTEP_SUBDIVISIONS, dt );
+    let ltaCircuit = new LTACircuit( ltaResistors, ltaBatteries, ltaCapacitors, ltaInductors );
+    let circuitResult = ltaCircuit.solveWithSubdivisions( TIMESTEP_SUBDIVISIONS, dt );
 
-    ltaResistors.forEach( resistorAdapter => {
-      const circuitElement = resistorMap.get( resistorAdapter )!;
-      if ( circuitElement instanceof LightBulb && circuitElement.isReal ) {
+    if ( hasRealBulbs ) {
+      ltaResistors.forEach( resistorAdapter => {
+        const circuitElement = resistorMap.get( resistorAdapter )!;
+        if ( circuitElement instanceof LightBulb && circuitElement.isReal ) {
 
-        const logWithBase = ( value: number, base: number ) => Math.log( value ) / Math.log( base );
+          const logWithBase = ( value: number, base: number ) => Math.log( value ) / Math.log( base );
 
-        const dV = circuitResult.getFinalState().ltaSolution!.getVoltage( resistorAdapter.nodeId0, resistorAdapter.nodeId1 );
-        const V = Math.abs( dV );
+          const dV = circuitResult.getFinalState().ltaSolution!.getVoltage( resistorAdapter.nodeId0, resistorAdapter.nodeId1 );
+          const V = Math.abs( dV );
 
-        const base = 2;
+          const base = 2;
 
-        // I = ln(V)
-        // V=IR
-        // V=ln(V)R
-        // R = V/ln(V)
+          // I = ln(V)
+          // V=IR
+          // V=ln(V)R
+          // R = V/ln(V)
 
-        // Adjust so it looks good in comparison to a standard bulb
-        const coefficient = 3;
+          // Adjust so it looks good in comparison to a standard bulb
+          const coefficient = 3;
 
-        // shift by base so at V=0 the log is 1
-        resistorAdapter.resistance = LightBulb.REAL_BULB_COLD_RESISTANCE + coefficient * V / logWithBase( V + base, base );
-        circuitElement.resistanceProperty.value = resistorAdapter.resistance;
-      }
-    } );
+          // shift by base so at V=0 the log is 1
+          resistorAdapter.resistance = LightBulb.REAL_BULB_COLD_RESISTANCE + coefficient * V / logWithBase( V + base, base );
+          circuitElement.resistanceProperty.value = resistorAdapter.resistance;
+        }
+      } );
 
+      // Solve the system
+      ltaCircuit = new LTACircuit( ltaResistors, ltaBatteries, ltaCapacitors, ltaInductors );
+      circuitResult = ltaCircuit.solveWithSubdivisions( TIMESTEP_SUBDIVISIONS, dt );
+    }
+
     // Apply the solutions from the analysis back to the actual Circuit
     ltaBatteries.forEach( batteryAdapter => {
       const circuitElement = voltageSourceMap.get( batteryAdapter )!;

[matthew-blackman]

@samreid's patch above appears to fix the fluctuation issue entirely. @arouinfar and I will investigate to confirm that the real bulb's current-voltage curve matches the intended logarithmic behavior. If the IV curve looks good, we believe that this issue is nearing completion.

[matthew-blackman]

@samreid and I tested the above patch and completed the following tests:

1. We compared the behavior to production - it matches perfectly
2. We plotted I vs V for a real bulb to confirm logarithmic relationship

Since the updated behavior does not fluctuate and matches the desired behavior, we agree that this issue can be closed after documenting the updates to LinearTransitAnalysis.

[matthew-blackman]

@arouinfar this is looking good on my end. Feel free to close this issue if it looks good to you.

[arouinfar]

@samreid @matthew-blackman looks great! I no longer see instability in the real bulb's currentProperty and I confirmed we are maintaining a logarithmic relationship between I and V.