view sim.py @ 11:8d34a9eec184

Split out cost function so we can re-run it on DB data to allow faster adjustment.
author Daniel O'Connor <darius@dons.net.au>
date Sat, 18 Nov 2023 10:37:35 +1030
parents 2832aefd442c
children 49939e179c86
line wrap: on
line source

import datetime
import pathlib
from scipy.optimize import differential_evolution
from spicelib.simulators.ltspice_simulator import LTspice
from spicelib.log.ltsteps import LTSpiceLogReader
from spicelib.editor.spice_editor import SpiceEditor
#from spicelib.sim.sim_runner import SimRunner
#import shlex
import sqlite3
import subprocess

def simulate(src, simexe, simflags, rundir, runname, params):
    src = pathlib.Path(src)
    runname = pathlib.Path(runname)
    assert(runname.suffix == '.net')
    rundir = pathlib.Path(rundir)

    e = SpiceEditor(src)
    for k in params:
        e.set_parameter(k, params[k])

    #s = SimRunner(simulator = LTspice, output_folder = 'tmp', verbose = True)
    #s.simulator.spice_exe = [simexe]
    #raw_file, log_file = s.run_now(e, run_filename = runname)

    e.write_netlist(rundir / runname)

    cmd = [simexe, '-b']
    cmd.extend(simflags)
    cmd.append(runname.name)
    #print(' '.join(map(shlex.quote, cmd)))
    then = datetime.datetime.now()
    p = subprocess.Popen(cmd, cwd = rundir)
    p.communicate()
    now = datetime.datetime.now()
    taken = now - then

    if p.returncode != 0:
        raise Exception('failed')

    #rawpath = rundir / runname.with_suffix('.raw')
    #raw = spicelib.RawRead(rawpath)
    logpath = rundir / runname.with_suffix('.log')
    log = LTSpiceLogReader(logpath)

    power = log.get_measure_value('pout')
    eff = log.get_measure_value('efficiency')
    thd = log.fourier['V(rfout)'][0].thd
    ipeak_u2 = log.get_measure_value('ipeak_u2')
    ipeak_u5 = log.get_measure_value('ipeak_u5')
    ipeak = max(ipeak_u2, ipeak_u5)

    return then, taken, power, eff, thd, ipeak

def calccost(power, eff, thd, ipeak):
    # Calculate the cost
    tpwr = 1500
    tthd = 2
    teff = 90
    imax = 11

    cost = 0
    if power < tpwr * 0.80:
        cost += 100
    elif power < tpwr:
        cost += (tpwr - power) / tpwr / 100

    if thd > 5 * tthd:
        cost += 100
    else:
        thdinv = 100 - thd
        tthdinv = 100 - tthd
        if thdinv < tthdinv:
            cost += (tthdinv - thdinv) / tthdinv
    if eff < teff:
        cost += (teff - eff) / teff

    if ipeak > imax:
        cost += (ipeak - imax) * 5

    return cost

def fn(v, cur, simexe, simflags, rundir, circ, ctr):
    '''Called by differential_evolution, v contains the evolved parameters, the rest are passed in from the args parameter'.
    Returns the cost value which is being minimised'''

    # Get parameters for run
    duty, c1, c2, l1, l2 = v

    # Check if this combination has already been tried
    cur.execute('SELECT run, power, efficiency, thd, ipeak FROM GAN190 WHERE duty = ? AND c1 = ? AND c2 = ? AND l1 = ? AND l2 = ?', (duty, c1, c2, l1, l2))
    tmp = cur.fetchone()
    if tmp is not None:
        run, power, eff, thd, ipeak = tmp
        # Recalculate the cost since it is cheap and might have been tweaked since the original run
        cost = calccost(power, eff, thd, ipeak)
        print(f'Found run {run:3d}: Duty {duty:3.0f}%, C1 {c1:3.0f}pF, C2 {c2:3.0f}pF, L1 {l1:3.0f}uH, L2 {l2:4.0f}nH -> Power: {power:6.1f}W Efficiency: {eff:5.1f}% THD: {thd:5.1f}% IPeak: {ipeak:4.1f}A Cost: {cost:6.2f}')
        return cost

    # Get next run number
    run = next(ctr)

    # Run the simulation
    # Need to convert units to suit
    runname = pathlib.Path(circ).stem + '-' + str(run) + '.net'
    then, taken, power, eff, thd, ipeak = simulate(circ, simexe, simflags, rundir, runname,
                                                   {'dutypct' : duty, 'C1' : c1 * 1e-12, 'C2' : c2 * 1e-12, 'L1' : l1 * 1e-6, 'L2' : l2 * 1e-9})
    # Calculate the cost
    cost = calccost(power, eff, thd, ipeak)

    # Log & save the results
    print(f'Run {run:3d}: Duty {duty:3.0f}%, C1 {c1:3.0f}pF, C2 {c2:3.0f}pF, L1 {l1:3.0f}uH, L2 {l2:4.0f}nH -> Power: {power:6.1f}W Efficiency: {eff:5.1f}% THD: {thd:5.1f}% IPeak: {ipeak:4.1f}A Cost: {cost:6.2f}')
    taken = taken.seconds + taken.microseconds / 1e6
    cur.execute('INSERT INTO GAN190 (name, run, start, duration, duty, c1, c2, l1, l2, power, efficiency, thd, ipeak, cost) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)',
                 (runname, run, then, taken, duty, c1, c2, l1, l2, power, eff, thd, ipeak, cost))
    cur.execute('COMMIT')

    return cost

def ev():
    # Bounds for parameters
    bounds = [(10, 80),
              (1, 20),
              (10, 300),
              (1, 20),
              (10, 500)]
    # Initial solution
    x0 = [36, 10, 155, 5, 140]
    dbh = sqlite3.connect('results.db')
    cur = dbh.cursor()
    cur.execute('''
CREATE TABLE IF NOT EXISTS GAN190 (
    name	TEXT,		-- Circuit name
    run		INTEGER,	-- Run number
    start	DATETIME,	-- Datetime run started
    duration	REAL,		-- Length of run (seconds)
    duty	REAL,		-- Duty cyle (%)
    c1		REAL,		-- Value of C1 (pF)
    c2		REAL,		-- Value of C2 (pF)
    l1		REAL,		-- Value of L1 (uH)
    l2		REAL,		-- Value of L2 (nH)
    power	REAL,		-- Measured power (Watts)
    efficiency	REAL,		-- Measured efficiency (%)
    thd		REAL,		-- Total harmonic distortion (%)
    ipeak	REAL,		-- Peak drain current (A)
    cost	REAL		-- Calculated cost metric
);''')
    def ctr():
        i = 0
        while True:
            yield i
            i += 1
    return differential_evolution(fn, bounds, x0 = x0,
                                  args = (cur, '/Users/oconnd1/bin/runltspice', ['-alt'], 'tmp', 'pa-GAN190-PP-nodriver.net', ctr()),
                                  integrality = True, disp = True, seed = 12345)