diff --git a/docs/notebooks/neuralnet/mnist_example_convolutional.ipynb b/docs/notebooks/neuralnet/mnist_example_convolutional.ipynb
index 8365ade9..e4e5fc02 100644
--- a/docs/notebooks/neuralnet/mnist_example_convolutional.ipynb
+++ b/docs/notebooks/neuralnet/mnist_example_convolutional.ipynb
@@ -53,7 +53,7 @@
     "#omlt for interfacing our neural network with pyomo\n",
     "from omlt import OmltBlock\n",
     "from omlt.neuralnet import FullSpaceNNFormulation\n",
-    "from omlt.io.onnx import write_onnx_model_with_bounds, load_onnx_neural_network_with_bounds"
+    "from omlt.io.onnx_reader import write_onnx_model_with_bounds, load_onnx_neural_network_with_bounds"
    ]
   },
   {
@@ -659,4 +659,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 2
-}
+}
\ No newline at end of file
diff --git a/docs/notebooks/neuralnet/mnist_example_dense.ipynb b/docs/notebooks/neuralnet/mnist_example_dense.ipynb
index 052059bb..e7f5b79d 100644
--- a/docs/notebooks/neuralnet/mnist_example_dense.ipynb
+++ b/docs/notebooks/neuralnet/mnist_example_dense.ipynb
@@ -52,7 +52,7 @@
     "#omlt for interfacing our neural network with pyomo\n",
     "from omlt import OmltBlock\n",
     "from omlt.neuralnet import FullSpaceNNFormulation\n",
-    "from omlt.io.onnx import write_onnx_model_with_bounds, load_onnx_neural_network_with_bounds"
+    "from omlt.io.onnx_reader import write_onnx_model_with_bounds, load_onnx_neural_network_with_bounds"
    ]
   },
   {
@@ -742,4 +742,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 2
-}
+}
\ No newline at end of file
diff --git a/src/omlt/io/__init__.py b/src/omlt/io/__init__.py
index 0b45abf7..67a5feac 100644
--- a/src/omlt/io/__init__.py
+++ b/src/omlt/io/__init__.py
@@ -1,2 +1,2 @@
-from omlt.io.onnx import load_onnx_neural_network, write_onnx_model_with_bounds, load_onnx_neural_network_with_bounds
+from omlt.io.onnx_reader import load_onnx_neural_network, write_onnx_model_with_bounds, load_onnx_neural_network_with_bounds
 from omlt.io.keras_reader import load_keras_sequential
diff --git a/src/omlt/io/onnx_parser.py b/src/omlt/io/onnx_parser.py
index 34de5ec1..6ac7a9d1 100644
--- a/src/omlt/io/onnx_parser.py
+++ b/src/omlt/io/onnx_parser.py
@@ -10,7 +10,7 @@
 )
 
 
-_ACTIVATION_OP_TYPES = ["Relu", "Sigmoid", "LogSoftmax"]
+_ACTIVATION_OP_TYPES = ["Relu", "Sigmoid", "LogSoftmax", "Tanh"]
 
 
 class NetworkParser:
@@ -181,6 +181,10 @@ def _consume_dense_nodes(self, node, next_nodes):
             node_biases = self._initializers[in_1]
 
         assert len(node_weights.shape) == 2
+
+        # Flatten biases array as some APIs (scikit) store biases as (1, n) instead of (n,)
+        node_biases = node_biases.flatten()
+
         assert node_weights.shape[1] == node_biases.shape[0]
         assert len(node.output) == 1
 
@@ -339,7 +343,12 @@ def _consume_reshape_nodes(self, node, next_nodes):
         assert len(node.input) == 2
         [in_0, in_1] = list(node.input)
         input_layer = self._node_map[in_0]
-        new_shape = self._constants[in_1]
+
+        if in_1 in self._constants:
+            new_shape = self._constants[in_1]
+        else:
+            new_shape = self._initializers[in_1]
+
         output_size = np.empty(input_layer.output_size).reshape(new_shape).shape
         transformer = IndexMapper(input_layer.output_size, list(output_size))
         self._node_map[node.output[0]] = (transformer, input_layer)
diff --git a/src/omlt/io/onnx.py b/src/omlt/io/onnx_reader.py
similarity index 100%
rename from src/omlt/io/onnx.py
rename to src/omlt/io/onnx_reader.py
diff --git a/src/omlt/io/sklearn_reader.py b/src/omlt/io/sklearn_reader.py
new file mode 100644
index 00000000..717a08fa
--- /dev/null
+++ b/src/omlt/io/sklearn_reader.py
@@ -0,0 +1,67 @@
+from skl2onnx.common.data_types import FloatTensorType
+from skl2onnx import convert_sklearn
+from sklearn.preprocessing import MinMaxScaler, MaxAbsScaler, StandardScaler
+from sklearn.preprocessing import RobustScaler
+from omlt.scaling import OffsetScaling
+from omlt.io.onnx_reader import load_onnx_neural_network
+import onnx
+
+def parse_sklearn_scaler(sklearn_scaler):
+
+    if isinstance(sklearn_scaler, StandardScaler):
+        offset = sklearn_scaler.mean_
+        factor = sklearn_scaler.scale_
+
+    elif isinstance(sklearn_scaler, MaxAbsScaler):
+        factor = sklearn_scaler.scale_
+        offset = factor*0
+
+    elif isinstance(sklearn_scaler, MinMaxScaler):
+        factor = sklearn_scaler.data_max_ - sklearn_scaler.data_min_
+        offset = sklearn_scaler.data_min_
+
+    elif isinstance(sklearn_scaler, RobustScaler):
+        factor = sklearn_scaler.scale_
+        offset = sklearn_scaler.center_
+
+    else:
+        raise(ValueError("Scaling object provided is not currently supported. Only linear scalers are supported."
+                         "Supported scalers include StandardScaler, MinMaxScaler, MaxAbsScaler, and RobustScaler"))
+
+    return offset, factor
+
+def convert_sklearn_scalers(sklearn_input_scaler, sklearn_output_scaler):
+
+    #Todo: support only scaling input or output?
+
+    offset_inputs, factor_inputs = parse_sklearn_scaler(sklearn_input_scaler)
+    offset_outputs, factor_ouputs = parse_sklearn_scaler(sklearn_output_scaler)
+
+    return OffsetScaling(offset_inputs=offset_inputs, factor_inputs=factor_inputs,
+                         offset_outputs=offset_outputs, factor_outputs=factor_ouputs)
+
+def load_sklearn_MLP(model, scaling_object=None, input_bounds=None, initial_types=None):
+
+    # Assume float inputs if no types are supplied to the model
+    if initial_types is None:
+        initial_types = [('float_input', FloatTensorType([None, model.n_features_in_]))]
+
+    onx = convert_sklearn(model, initial_types=initial_types, target_opset=12)
+
+    # Remove initial cast layer created by sklearn2onnx
+    graph = onx.graph
+    node1 = graph.node[0]
+    graph.node.remove(node1)
+    new_node = onnx.helper.make_node(
+        'MatMul',
+        name="MatMul",
+        inputs=['float_input', 'coefficient'],
+        outputs=['mul_result']
+    )
+    graph.node.insert(0, new_node)
+
+    # Replace old MatMul node with new node with correct input name
+    node2 = graph.node[1]
+    graph.node.remove(node2)
+
+    return load_onnx_neural_network(onx, scaling_object, input_bounds)
diff --git a/tests/io/test_onnx_parser.py b/tests/io/test_onnx_parser.py
index 402c7e39..257a0186 100644
--- a/tests/io/test_onnx_parser.py
+++ b/tests/io/test_onnx_parser.py
@@ -2,7 +2,7 @@
 import onnx
 import numpy as np
 
-from omlt.io.onnx import load_onnx_neural_network
+from omlt.io.onnx_reader import load_onnx_neural_network
 
 
 def test_linear_131(datadir):
diff --git a/tests/io/test_sklearn.py b/tests/io/test_sklearn.py
new file mode 100644
index 00000000..3e6f23dd
--- /dev/null
+++ b/tests/io/test_sklearn.py
@@ -0,0 +1,164 @@
+from omlt.scaling import OffsetScaling
+from omlt.block import OmltBlock
+from omlt.io.sklearn_reader import convert_sklearn_scalers
+from omlt.io.sklearn_reader import load_sklearn_MLP
+from omlt.neuralnet import FullSpaceNNFormulation
+from sklearn.preprocessing import MinMaxScaler, MaxAbsScaler, StandardScaler, RobustScaler
+from scipy.stats import iqr
+from pyomo.environ import *
+import numpy as np
+import json
+import pickle
+
+def test_sklearn_scaler_conversion():
+    X = np.array(
+        [[42, 10, 29],
+         [12, 19, 15]]
+    )
+
+    Y = np.array(
+        [[1, 2],
+         [3, 4]]
+    )
+
+    # Create sklearn scalers
+    xMinMax = MinMaxScaler()
+    xMaxAbs = MaxAbsScaler()
+    xStandard = StandardScaler()
+    xRobust = RobustScaler()
+
+    yMinMax = MinMaxScaler()
+    yMaxAbs = MaxAbsScaler()
+    yStandard = StandardScaler()
+    yRobust = RobustScaler()
+
+    sklearn_scalers = [(xMinMax, yMinMax), (xMaxAbs, yMaxAbs), (xStandard, yStandard), (xRobust, yRobust)]
+    for scalers in sklearn_scalers:
+        scalers[0].fit(X)
+        scalers[1].fit(Y)
+
+    # Create OMLT scalers using OMLT function
+    MinMaxOMLT = convert_sklearn_scalers(xMinMax, yMinMax)
+    MaxAbsOMLT = convert_sklearn_scalers(xMaxAbs, yMaxAbs)
+    StandardOMLT = convert_sklearn_scalers(xStandard, yStandard)
+    RobustOMLT = convert_sklearn_scalers(xRobust, yRobust)
+
+    omlt_scalers = [MinMaxOMLT, MaxAbsOMLT, StandardOMLT, RobustOMLT]
+
+    # Generate test data
+    x = {0: 10, 1: 29, 2: 42}
+    y = {0: 2, 1: 1}
+
+    # Test Scalers
+    for i in range(len(omlt_scalers)):
+        x_s_omlt = omlt_scalers[i].get_scaled_input_expressions(x)
+        y_s_omlt = omlt_scalers[i].get_scaled_output_expressions(y)
+
+        x_s_sklearn = sklearn_scalers[i][0].transform([list(x.values())])[0]
+        y_s_sklearn = sklearn_scalers[i][1].transform([list(y.values())])[0]
+
+        np.testing.assert_almost_equal(list(x_s_omlt.values()), list(x_s_sklearn))
+        np.testing.assert_almost_equal(list(y_s_omlt.values()), list(y_s_sklearn))
+
+def test_sklearn_offset_equivalence():
+    X = np.array(
+        [[42, 10, 29],
+         [12, 19, 15]]
+    )
+
+    Y = np.array(
+        [[1, 2],
+         [3, 4]]
+    )
+
+    # Get scaling factors for OffsetScaler
+    xmean = X.mean(axis=0)
+    xstd = X.std(axis=0)
+    xmax = X.max(axis=0)
+    absxmax = abs(X).max(axis=0)
+    xmin = X.min(axis=0)
+    xminmax = xmax-xmin
+    xmedian = np.median(X, axis=0)
+    xiqr = iqr(X, axis=0)
+
+    ymean = Y.mean(axis=0)
+    ystd = Y.std(axis=0)
+    ymax = Y.max(axis=0)
+    absymax = abs(Y).max(axis=0)
+    ymin = Y.min(axis=0)
+    yminmax = ymax-ymin
+    ymedian = np.median(Y, axis=0)
+    yiqr = iqr(Y, axis=0)
+
+    # Create sklearn scalers
+    xMinMax = MinMaxScaler()
+    xMaxAbs = MaxAbsScaler()
+    xStandard = StandardScaler()
+    xRobust = RobustScaler()
+
+    yMinMax = MinMaxScaler()
+    yMaxAbs = MaxAbsScaler()
+    yStandard = StandardScaler()
+    yRobust = RobustScaler()
+
+    sklearn_scalers = [(xMinMax, yMinMax), (xMaxAbs, yMaxAbs), (xStandard, yStandard), (xRobust, yRobust)]
+    for scalers in sklearn_scalers:
+        scalers[0].fit(X)
+        scalers[1].fit(Y)
+
+    # Create OMLT scalers manually
+    MinMaxOMLT = OffsetScaling(offset_inputs=xmin, factor_inputs=xminmax, offset_outputs=ymin, factor_outputs=yminmax)
+    MaxAbsOMLT = OffsetScaling(offset_inputs=[0]*3, factor_inputs=absxmax, offset_outputs=[0]*2, factor_outputs=absymax)
+    StandardOMLT = OffsetScaling(offset_inputs=xmean, factor_inputs=xstd, offset_outputs=ymean, factor_outputs=ystd)
+    RobustOMLT = OffsetScaling(offset_inputs=xmedian, factor_inputs=xiqr, offset_outputs=ymedian, factor_outputs=yiqr)
+
+    omlt_scalers = [MinMaxOMLT, MaxAbsOMLT, StandardOMLT, RobustOMLT]
+
+    # Generate test data
+    x = {0: 10, 1: 29, 2: 42}
+    y = {0: 2, 1: 1}
+
+    # Test Scalers
+    for i in range(len(omlt_scalers)):
+        x_s_omlt = omlt_scalers[i].get_scaled_input_expressions(x)
+        y_s_omlt = omlt_scalers[i].get_scaled_output_expressions(y)
+
+        x_s_sklearn = sklearn_scalers[i][0].transform([list(x.values())])[0]
+        y_s_sklearn = sklearn_scalers[i][1].transform([list(y.values())])[0]
+
+        np.testing.assert_almost_equal(list(x_s_omlt.values()), list(x_s_sklearn))
+        np.testing.assert_almost_equal(list(y_s_omlt.values()), list(y_s_sklearn))
+
+def test_sklearn_model(datadir):
+    nn_names = ["sklearn_identity_131", "sklearn_logistic_131", "sklearn_tanh_131"]
+
+    # Test each nn
+    for nn_name in nn_names:
+        nn = pickle.load(open(datadir.file(nn_name+".pkl"), 'rb'))
+
+        with open(datadir.file(nn_name+"_bounds"), 'r') as f:
+            bounds = json.load(f)
+
+            # Convert to omlt format
+            xbounds = {int(i): tuple(bounds[i]) for i in bounds}
+
+        net = load_sklearn_MLP(nn, input_bounds=xbounds)
+        formulation = FullSpaceNNFormulation(net)
+
+        model = ConcreteModel()
+        model.nn = OmltBlock()
+        model.nn.build_formulation(formulation)
+
+        @model.Objective()
+        def obj(mdl):
+            return 1
+
+        x = [(xbounds[i][0]+xbounds[i][1])/2.0 for i in range(2)]
+        for i in range(len(x)):
+            model.nn.inputs[i].fix(x[i])
+
+        result = SolverFactory("ipopt").solve(model, tee=False)
+        yomlt = [value(model.nn.outputs[0]), value(model.nn.outputs[1])]
+
+        ysklearn = nn.predict([x])[0]
+        np.testing.assert_almost_equal(list(yomlt), list(ysklearn))
\ No newline at end of file
diff --git a/tests/models/sklearn_identity_131.pkl b/tests/models/sklearn_identity_131.pkl
new file mode 100644
index 00000000..f5889b2d
Binary files /dev/null and b/tests/models/sklearn_identity_131.pkl differ
diff --git a/tests/models/sklearn_identity_131_bounds b/tests/models/sklearn_identity_131_bounds
new file mode 100644
index 00000000..007637c7
--- /dev/null
+++ b/tests/models/sklearn_identity_131_bounds
@@ -0,0 +1 @@
+{"0": [-3.2412673400690726, 2.3146585666735087], "1": [-1.9875689146008928, 3.852731490654721]}
\ No newline at end of file
diff --git a/tests/models/sklearn_logistic_131.pkl b/tests/models/sklearn_logistic_131.pkl
new file mode 100644
index 00000000..9cb2b3d2
Binary files /dev/null and b/tests/models/sklearn_logistic_131.pkl differ
diff --git a/tests/models/sklearn_logistic_131_bounds b/tests/models/sklearn_logistic_131_bounds
new file mode 100644
index 00000000..007637c7
--- /dev/null
+++ b/tests/models/sklearn_logistic_131_bounds
@@ -0,0 +1 @@
+{"0": [-3.2412673400690726, 2.3146585666735087], "1": [-1.9875689146008928, 3.852731490654721]}
\ No newline at end of file
diff --git a/tests/models/sklearn_tanh_131.pkl b/tests/models/sklearn_tanh_131.pkl
new file mode 100644
index 00000000..b0732822
Binary files /dev/null and b/tests/models/sklearn_tanh_131.pkl differ
diff --git a/tests/models/sklearn_tanh_131_bounds b/tests/models/sklearn_tanh_131_bounds
new file mode 100644
index 00000000..007637c7
--- /dev/null
+++ b/tests/models/sklearn_tanh_131_bounds
@@ -0,0 +1 @@
+{"0": [-3.2412673400690726, 2.3146585666735087], "1": [-1.9875689146008928, 3.852731490654721]}
\ No newline at end of file
diff --git a/tests/neuralnet/test_onnx.py b/tests/neuralnet/test_onnx.py
index 2258b975..73e040d7 100644
--- a/tests/neuralnet/test_onnx.py
+++ b/tests/neuralnet/test_onnx.py
@@ -1,5 +1,5 @@
 import tempfile
-from omlt.io.onnx import load_onnx_neural_network, load_onnx_neural_network_with_bounds, write_onnx_model_with_bounds
+from omlt.io.onnx_reader import load_onnx_neural_network, load_onnx_neural_network_with_bounds, write_onnx_model_with_bounds
 import onnx
 import onnxruntime as ort
 import numpy as np
diff --git a/tests/neuralnet/test_relu.py b/tests/neuralnet/test_relu.py
index 354dd322..3f6e9669 100644
--- a/tests/neuralnet/test_relu.py
+++ b/tests/neuralnet/test_relu.py
@@ -2,7 +2,7 @@
 import numpy as np
 
 from omlt.block import OmltBlock
-from omlt.io.onnx import load_onnx_neural_network_with_bounds
+from omlt.io.onnx_reader import load_onnx_neural_network_with_bounds
 from omlt.neuralnet import FullSpaceNNFormulation, ReluBigMFormulation, ReluComplementarityFormulation, ReluPartitionFormulation
 from omlt.neuralnet.activations import ComplementarityReLUActivation